Modeling flow through non linear coarse porous media by a finite difference scheme, with application to rockfill dams overtopping

==Resumen==

La circulación del agua a través de medios porosos gruesos, como la escollera, la grava o incluso la arena de mayor tamaño, tiene la particularidad de exhibir una relación no lineal entre el gradiente hidráulico y la velocidad de filtración. Uno de los problemas de ingeniería en que se observa esto es cuando se produce un sobrevertido en una presa de escollera y el agua, pasando sobre su coronación, se infiltra en el espaldón y circula por su interior. La modelación numérica del problema se resuelve aquí mediante un esquema de diferencias finitas, que se implementa en el código MNLEE.RDS para su aplicación al sobrevertido en presas de escollera. Además, se proporciona un procedimiento para resolver el acoplamiento del movimiento del agua en lámina libre sobre el talud de la presa y de filtración por el interior del espaldón, mediante la introducción del concepto «ley de intercambio de caudales».

==Abstract==

Flow through coarse porous media, as rockfill, gravel or sand of big size, exhibit a non linear relation between hydraulic gradient and seepage velocity. This behavior is observed when water flows over the crest of a rockfill dam, seeps into the downstream shoulder and moves through the rockfill voids. A finite difference scheme is used here to solve this problem through numerical modeling. It is implemented in the code MNLEE.RDS for solving dam overtopping problems. In addition, a procedure is provided to solve the coupled problem of flow over and through the rockfill dam. For that purpose, the concept of «law of flow exchange»is introduced.

==Palabras clave==

Filtración ; Medios porosos ; Diferencias finitas ; Escollera ; Presa ; Sobrevertido

==Keywords==

Seepage ; Porous media ; Finite differences ; Rockfill ; Dam ; Overtopping

==1. Introducción==

Al realizar estudios de filtración del agua a través de medios porosos finos, como el núcleo de arcilla de una presa de materiales sueltos, se considera una relación lineal entre el gradiente de carga y la velocidad de filtración, de acuerdo con la conocida ley de Darcy. La experiencia ha demostrado que en medios granulares gruesos (escolleras, gravas y arenas gruesas) esta ley no es válida, pues dicha relación no es lineal. En este sentido podemos hablar de medios porosos no lineales. Como ejemplo de ello podemos citar la filtración en acuíferos de grano grueso o en materiales granulares utilizados como elemento de drenaje. También es el caso de la circulación de agua que se establece por el interior del espaldón de una presa de escollera cuando excepcionalmente, debido a la insuficiencia de capacidad de su aliviadero, durante una avenida el agua llega a pasar sobre su coronación. Actualmente, las distintas administraciones hidráulicas de gran número de países están abordando el estudio y mejora de la seguridad de sus presas, teniendo en este contexto una importancia singular el problema del sobrevertido en presas de materiales sueltos, pues ha sido causa de gran número de roturas en este tipo de presas. Los principales mecanismos involucrados en la rotura de una presa de escollera por sobrevertido son el deslizamiento en masa y el arrastre de partículas. Además, existe el peligro de que se produzca migración interna de partículas. Por tanto, para evaluar la seguridad o para diseñar una presa resistente a este fenómeno, es condición previa conocer las presiones intersticiales y los gradientes hidráulicos y, por tanto, resolver el problema de filtración en un medio granular no lineal. Para ello, es necesario conocer la fórmula de resistencia que expresa la antes referida relación gradiente/velocidad y, a partir de ella, establecer una ecuación de campo que caracterice el fenómeno de filtración en este tipo de medios granulares gruesos. Las condiciones de frontera particularizarán el problema a resolver, de modo que el problema quede adecuadamente planteado. La resolución por métodos numéricos implica la discretización del dominio de filtración y la resolución del sistema de ecuaciones resultante. Se expone en apartados sucesivos cada uno de los aspectos señalados. Debe llamarse la atención sobre la necesidad de que la relación entre el tamaño del dominio de filtración y el tamaño de partícula sea suficientemente alta para asegurar la aplicabilidad de los procedimientos diferenciales que aquí se emplean. Esta consideración, que siempre aplica, puede resultar condicionante en el caso de que el tamaño de grano del material poroso sea grande, como sucede en la escollera.

==2. Fórmula de resistencia==

Diversos autores han propuesto leyes alternativas para medios no lineales, siendo las más usadas las leyes cuadrática y exponencial.

Ley cuadrática:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i=a\quad V_d+b\quad V_d^2\quad a,b:\quad constantes.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

Ley exponencial:

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i=c\quad V_d^m\quad c,m:\quad constantes\quad caracter\overset{\textasciiacute}{\imath }sticas\quad del\quad material.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

Siendo ''i''  el gradiente hidráulico y ''V''<sub>''d''</sub>  la velocidad de filtración. Esta última se define como la relación entre el caudal filtrado y la sección de paso total, incluidos huecos y esqueleto sólido.      

George y Hansen (1992)[[#bib0005|[1]]]  obtuvieron unas fórmulas que permiten determinar los coeficientes de la fórmula cuadrática a partir del coeficiente y el exponente de la fórmula exponencial y viceversa.      

Conversión de <math display="inline">i=c\quad V_d^m</math>  a <math display="inline">i=a\quad V_d+b\quad V_d^2</math>

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a=\frac{12\quad c(V_{d,max})^{m-1}(2-m)}{\left(m+2)(m+3\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>b=\frac{20\quad c(V_{d,max})^{m-2}(m-1)}{\left(m+2)(m+3\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

Siendo ''i''  = 0 para ''V''<sub>''d''</sub>  = 0 y ''V''<sub>''d'' ,''max''</sub>  el valor máximo de la velocidad de filtración en el intervalo considerado.      

Conversión de <math display="inline">i=a\quad V_d+b\quad V_d^2</math>  a <math display="inline">i=c\quad V_d^m</math>

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>m=\frac{5\quad a+6\quad b\quad V_{d,max}}{5\quad a+3\quad b\quad V_{d,max}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>c=\frac{\left(5\quad a+4\quad b\quad V_{d,max})(4\quad a+3\quad b\quad V_{d,max}\right)}{4(5\quad a+3\quad b\quad V_{d,max})(V_{d,max})^{m-1}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

La modelación se ha realizado partiendo de una fórmula de resistencia de tipo exponencial. No obstante, si se dispone de una fórmula de tipo cuadrático para un material determinado, basta con la conversión mediante las anteriores fórmulas de George y Hansen. A partir de resultados obtenidos empíricamente, Parkin, Trollope y Lawson (1966)[[#bib0010|[2]]]  definieron un nomograma que permite obtener el coeficiente ''c''  de la fórmula exponencial a partir del índice de huecos y la superficie específica del medio granular considerado. Así mismo, llegaron a la conclusión de que el exponente de dicha fórmula varía poco en torno al valor 1,85 próximo al valor 2 correspondiente a movimiento totalmente turbulento, y alejado del valor 1 correspondiente a movimiento laminar. La causa de que el exponente de la fórmula exponencial tenga un valor próximo a 2, pero sin llegar al mismo, ha sido objeto de controversia. Con independencia de cuál sea la explicación, es un hecho comprobado que el valor de dicho exponente se encuentra, en efecto, en dicho rango de variación.      

==3. Movimiento del agua. Planteamiento general del problema==

Una vez conocida la fórmula de resistencia, puede plantearse el estudio del movimiento del agua a través de la escollera. El análisis diferencial conduce a un modelo matemático diferencial que se plantea como un problema de campo. El análisis se limitará al caso de movimiento permanente, bidimensional de un fluido viscoso e incompresible, el agua, en un medio poroso homogéneo e isótropo. En primer lugar es preciso definir el problema. Para ello debe conocerse la ecuación de campo que caracteriza el fenómeno físico: movimiento del agua a través de un medio poroso determinado, y las condiciones de frontera, que particularizan el problema para un caso concreto. Como resultado final del cálculo, se obtiene un campo de presiones y un campo de velocidades, si bien la ecuación de campo puede plantearse en términos de otras magnitudes, típicamente su potencial. En cualquier caso, una vez resuelto el problema, es posible conocer la presión y la velocidad en cualquier punto del dominio de filtración y, a través de la velocidad, el caudal de agua filtrado. Por tanto, el planteamiento del problema consiste en la definición de la ecuación de campo y de las condiciones de frontera. Cuando hablamos de resolver el problema, nos referimos a la determinación de las presiones y velocidades en cada punto del medio. Se dice que el problema está bien planteado cuando la solución existe y es única. Por otra parte, al plantear un problema de campo es necesario garantizar que la ecuación de campo se cumple en todo el dominio. El incumplimiento de esta condición invalida el planteamiento del problema.

==4. Ecuación de campo==

Consideremos una partícula de medio poroso de dimensiones suficientemente reducidas en relación con el dominio de filtración para que los procedimientos diferenciales sean aplicables: ''dxdy'' , y unos ejes de coordenadas ortogonales (''x'' , ''y'' ). Sea el vector velocidad de filtración en un punto:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>V=u\boldsymbol{\mbox{x}}+v\boldsymbol{\mbox{y}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

siendo '''x'''  e '''y'''  vectores unitarios de dirección la de los ejes de coordenadas ''x''  e ''y''  respectivamente y sentido el definido como positivo. ( [[#fig0005|fig. 1]] )

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr1.jpg|center|351px|Partícula diferencial del medio poroso.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 1.

Partícula diferencial del medio poroso.

</span>
|}

El caudal de agua que entra en la partícula de medio poroso es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>udy+vdx</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

El caudal que sale es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(u+\frac{\delta u}{\delta x}\right)dy+\left(v+\right. </math><math>\left. \frac{\delta v}{\delta y}\right)dx</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

Admitiendo que ni el agua ni los elementos sólidos sufren deformación volumétrica y que la partícula está saturada, el caudal de agua que entra en ella tiene que ser igual que el que sale (ley de conservación de la masa o principio de continuidad), de donde se obtiene:

<span id='eq0050'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\delta u}{\delta x}+\frac{\delta v}{\delta y}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

===4.1. Medios porosos lineales===

Si se cumple la ley de Darcy, tendremos que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>V=Ki</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

siendo ''K''  el coeficiente de permeabilidad e '''i'''  el gradiente hidráulico. Por tanto:

<span id='eq0060'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u=K\frac{dh}{dx}\quad ;\quad v=K\frac{dh}{dy}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>h=z+\frac{p}{\gamma }\quad :carga\quad piezom\acute{e}trica</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Donde: ''z''  es la altura sobre el plano de referencia y ''p''  es la presión. En realidad, ''h''  es la carga hidráulica total, pero la carga de velocidad tiene en la práctica un valor muy bajo, despreciable frente al de la carga piezométrica. Sustituyendo  [[#eq0060|(12)]]  en [[#eq0050|(10)]]  tenemos que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\delta }{\delta x}\left(K\frac{\delta h}{\delta x}\right)+</math><math>\frac{\delta }{\delta y}\left(K\frac{\delta h}{\delta y}\right)=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

En un medio homogéneo e isótropo el coeficiente de permeabilidad es constante en todo el dominio. Por tanto:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\delta }^2h}{\delta x^2}+\frac{{\delta }^2h}{\delta y^2}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

Esta es la conocida ecuación de Laplace, que gobierna el movimiento del agua en un medio poroso cuando se cumple la ley de Darcy. Se trata, por tanto, de la ecuación de campo para movimiento lineal.

Si consideramos el operador laplaciano:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Delta =\frac{{\delta }^2}{\delta x^2}+\frac{{\delta }^2}{\delta y^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

la ecuación de Laplace puede expresarse:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Delta h=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

===4.2. Medios porosos no lineales===

Diversos autores han desarrollado ecuaciones de campo para medios porosos no lineales basándose en una fórmula de resistencia de tipo exponencial o cuadrático.

Engelund (1953)[[#bib0015|[3]]]  y otros, considerando una fórmula de resistencia cuadrática con término lineal nulo, han llegado a la siguiente ecuación de campo.

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\delta }^2h}{\delta V^2}+\frac{2}{V^2}\frac{{\delta }^2h}{\delta {\alpha }^2}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)
|}

siendo: ''h''  carga piezométrica, ''V''  la velocidad de filtración y ''α''  la inclinación del vector velocidad. A continuación veremos el desarrollo completo de las ecuaciones de campo de Parkin y de Volker para movimiento no lineal.      

====4.2.1. Ecuación de campo de Parkin====

La ecuación de campo de Parkin (1971)[[#bib0020|[4]]]  se basa en una ley exponencial:

<span id='eq0095'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i=c\quad V^m</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

admitiendo que los coeficientes ''c''  y ''m''  son constantes en el dominio de filtración. Estos coeficientes no varían de forma significativa en el rango de velocidades que se dan en la práctica, por lo que la hipótesis resulta aceptable. De  [[#eq0095|(18)]]  se obtiene que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>V=(i/c)^{1/m}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

Tanto la velocidad como la pérdida unitaria de energía se refieren a la dirección de avance del agua en el punto considerado (línea de corriente del movimiento en el medio poroso). En general, la relación puede expresarse en forma vectorial. Teniendo en cuenta que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{i}}=\boldsymbol{\mbox{grad}}\quad \quad H</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

donde, despreciando la carga de velocidad queda:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{V}}={\left(\frac{1}{c}\right)}^{1/m}{\left|\boldsymbol{\mbox{grad}}\quad \quad h\right|}^{\frac{1}{m}-1}\boldsymbol{\mbox{grad}}\quad \quad h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

Admitiendo que el coeficiente ''c''  es constante, puede expresarse:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{V}}={\left(\frac{1}{c}\right)}^{\frac{1}{m}-1}{\left|\boldsymbol{\mbox{grad}}\quad \quad h\right|}^{\frac{1}{m}-1}\boldsymbol{\mbox{grad}}\quad \quad \left(\frac{h}{c}\right)=</math><math>{\left|\boldsymbol{\mbox{grad}}\quad \quad \left(\frac{h}{c}\right)\right|}^{\frac{1}{m}-1}\boldsymbol{\mbox{grad}}\quad \quad \left(\frac{h}{c}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)
|}

y llamando:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{h}{c}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)
|}

Se tiene:

<span id='eq0125'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{V}}={\left|\boldsymbol{\mbox{grad}}\quad \quad \phi \right|}^{\frac{1}{m}-1}\boldsymbol{\mbox{grad}}\quad \quad \phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)
|}

Por definición del operador gradiente:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{grad}}\quad \quad \phi =\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+</math><math>\frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)
|}

Siendo '''x'''  e '''y'''  los vectores unitarios en las direcciones de los ejes cartesianos ''x''  e ''y''  respectivamente.      

Sustituyendo en [[#eq0125|(24)]] :

<span id='eq0135'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{V}}={\left|\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+\frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}\right|}^{\frac{1}{m}-1}\left(\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+\right. </math><math>\left. \frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)
|}

Por otra parte, tenemos que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\left|\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+\frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}\right|}^{\frac{1}{m}-1}=</math><math>{\left|\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+\frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}\right|}^{\frac{1}{2}\left(\frac{1}{m}-1\right)2}=</math><math>{\left({\left|\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+\frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}\right|}^2\right)}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}=</math><math>{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)
|}

Por tanto, sustituyendo en [[#eq0135|(26)]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{V}}={\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}\left(\frac{\delta \phi }{\delta x}\boldsymbol{\mbox{x}}+\right. </math><math>\left. \frac{\delta \phi }{\delta y}\boldsymbol{\mbox{y}}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)
|}

En consecuencia, las componentes del vector velocidad

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{V}}=u\boldsymbol{\mbox{x}}+v\boldsymbol{\mbox{y}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)
|}

son:

<span id='eq0155'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u=\frac{\delta \phi }{\delta x}{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{1/2(1/m-1)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)
|}

<span id='eq0160'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>v=\frac{\delta \phi }{\delta y}{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{1/2(1/m-1)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)
|}

Como hemos visto en el apartado anterior, la ecuación de continuidad aplicada a una partícula de medio poroso indeformable conduce a:

<span id='eq0165'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\delta u}{\delta x}+\frac{\delta v}{\delta y}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)
|}

Derivando las expresiones [[#eq0155|(30)]]  y [[#eq0160|(31)]]  respecto a ''x''  e ''y''  respectivamente, obtenemos:

<span id='eq0170'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\delta u}{\delta x}=\frac{{\delta }^2\phi }{\delta x^2}{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}+</math><math>\frac{\delta \phi }{\delta x}\frac{1}{2}\left(\frac{1}{m}-\right. </math><math>\left. 1\right).{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}.\left(2\frac{\delta \phi }{\delta x}\frac{{\delta }^2\phi }{\delta x^2}+\right. </math><math>\left. 2\frac{\delta \phi }{\delta y}\frac{{\delta }^2\phi }{\delta y\delta x}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)
|}

<span id='eq0175'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\delta v}{\delta y}=\frac{{\delta }^2\phi }{\delta y^2}{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}+</math><math>\frac{\delta \phi }{\delta y}\frac{1}{2}\left(\frac{1}{m}-\right. </math><math>\left. 1\right).{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}\left(2\frac{\delta \phi }{\delta x}\frac{{\delta }^2\phi }{\delta x\delta y}+\right. </math><math>\left. 2\frac{\delta \phi }{\delta y}\frac{{\delta }^2\phi }{\delta y^2}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)
|}

Sustituyendo las expresiones [[#eq0170|(33)]]  y [[#eq0175|(34)]]  en la ecuación de continuidad [[#eq0165|(32)]] , obtenemos la ecuación de campo de Parkin.

<span id='eq0180'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{{\delta }^2\phi }{\delta x^2}+\frac{{\delta }^2\phi }{\delta y^2}\right)\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+\right. </math><math>\left. {\left(\frac{\delta \phi }{\delta y}\right)}^2\right]+</math><math>\left(\frac{1}{m}-1\right)[{\left(\frac{\delta \phi }{\delta x}\right)}^2\frac{{\delta }^2\phi }{\delta x^2}+</math><math>2\frac{\delta \phi }{\delta x}\frac{\delta \phi }{\delta y}\frac{{\delta }^2\phi }{\delta x\delta y}+</math><math>{\left(\frac{\delta \phi }{\delta y}\right)}^2\frac{{\delta }^2\phi }{\delta y^2}]=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)
|}

donde recordemos que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{h}{c}\quad ;\quad h=z+\frac{p}{\gamma }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Parkin, Trollope y Lawson (1966)[[#bib0010|[2]]]  hacen notar que, supuesto constante el coeficiente ''c'' , la anterior ecuación de campo es equivalente a otra en la que ''ϕ''  se sustituya por ''h'' , para obtener directamente la carga piezométrica. En consecuencia, la carga piezométrica (''h'' ), y por tanto la presión (''p'' ), dependen únicamente del exponente (''m'' ) de la fórmula de resistencia. La ecuación de campo  [[#eq0180|(35)]]  se reduce a la ecuación de Laplace, válida para movimiento lineal, cuando ''m''  =1, o lo que es lo mismo, cuando se cumple la ley de Darcy.      

==5. Condiciones de frontera==

El fenómeno de filtración del agua a través del espaldón de una presa de escollera sometida a sobrevertido da origen a un problema de frontera mixto con superficie libre. Consideremos en primer lugar el movimiento del agua a través de una presa homogénea ([[#fig0010|fig. 2]] ). Las condiciones de frontera son:
* Tramo AB: la carga piezométrica es constante e igual a la altura de agua sobre la base (''H'' ). Se trata de una condición de tipo Dirichlet.                          
* Tramo DE: análogamente, la carga piezométrica es en este caso ''H'' ′. También es una condición de tipo Dirichlet.                          
* Tramo AE: admitiendo que la base es impermeable, el caudal que circula a través de este tramo de frontera es nulo. Por tanto, se trata de una condición tipo Neumann.
* Tramo BC: se trata de una superficie libre, cuya situación se desconoce a priori. La presión en todos sus puntos debe ser nula (referida a la presión atmosférica) y no puede circular agua a través de la misma. Se trata de una condición tipo Cauchy.
* Tramo CD: la carga piezométrica viene dada por la altura sobre el plano de referencia, ya que la presión es nula por encontrarse en contacto con la atmósfera. La condición es de tipo Dirichlet.  Si el agua supera la cota de coronación se producirá el sobrevertido. Cuando la altura de lámina vertiente (''H''<sub>''L''</sub> ) sea pequeña, la circulación superficial o exterior, sobre la coronación y el talud de aguas abajo, se reducirá a unos pocos metros. A partir de un determinado punto solo se producirá circulación interna a través del cuerpo de la presa ([[#fig0015|fig. 3]] a). A medida que aumenta la altura de lámina vertiente, el punto extremo de infiltración se desplazará hacia aguas abajo sobre el talud de la presa y el ''D''  lo hará hacia aguas arriba. Cuando la lámina sea suficientemente elevada, por ser el caudal de llegada al embalse netamente superior al que puede circular por el interior de la presa, se producirá movimiento superficial sobre toda la longitud del talud de aguas abajo ( [[#fig0015|fig. 3]] b). La condición de frontera en la coronación y el paramento de aguas abajo puede plantearse como de potencial impuesto, pero desconocido a priori, que es preciso determinar mediante tanteos sucesivos de modo semejante a como se procede con los problemas de superficie libre. En realidad, se trata de determinar la geometría de la superficie libre del agua que circula sobre el talud de la presa. Por tanto, para determinar la situación de una superficie libre se sigue el siguiente proceso:
* 1.° Se estima la situación de la superficie libre.
* 2.° Se calcula la solución al problema, obteniéndose la presión en la superficie libre.
* 3.° Si la presión no es nula en todos los puntos de la misma, se modifica su geometría, elevándola en la zona donde es positiva y reduciendo su cota donde es negativa.
* 4.° Se vuelve a realizar el cálculo y el proceso se repite hasta que se consigue que la presión sea nula en todos los puntos de la superficie libre.  En el caso de que exista circulación sobre el talud de la presa, la condición que debe cumplirse es la que se explica a continuación. Se consideran dos secciones verticales separadas una cierta distancia ''Δx'' , siendo:
* ''q''<sub>1(2),''e'' (''i'' )                                        </sub> : caudal que atraviesa la sección 1 [[#eq0010|(2)]]  circulando exteriormente (interiormente) sobre el talud de la presa (a través del cuerpo de la presa).                          
* ''Δq''<sub>1,2</sub> : caudal que pasa del exterior al interior de la presa (o viceversa, si el signo es negativo) entre las secciones 1 y 2.                                      Aplicando el principio de continuidad al recinto de medio poroso RSVU ([[#fig0020|fig. 4]] ), tenemos:

<span id='eq0190'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>q_{2,i}=q_{1,i}+\Delta q_{1,2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)
|}

El movimiento del agua sobre el talud es de caudal variable, dado que existe un intercambio de caudales con el interior de la presa.

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr2.jpg|center|375px|Presa homogénea sin elemento impermeable.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 2.

Presa homogénea sin elemento impermeable.

</span>
|}

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr3.jpg|center|360px|Fases de movimiento del agua sobre la presa a) Fase inicial b) Fase final.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 3.

Fases de movimiento del agua sobre la presa a) Fase inicial b) Fase final.

</span>
|}

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr4.jpg|center|334px|Caudales circulantes a través de las secciones paralelas.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 4.

Caudales circulantes a través de las secciones paralelas.

</span>
|}

El movimiento del agua sobre el talud y a través del cuerpo de la presa son fenómenos acoplados, no siendo posible en pura teoría calcular uno de ellos con independencia del otro. En efecto, para resolver el movimiento del agua a través de la presa es preciso conocer la carga piezométrica en el talud de aguas abajo, y por tanto el perfil del agua que circula sobre el mismo. Pero para determinar este perfil tenemos que conocer el caudal que se añade o sustrae en cada punto del talud, que se determina mediante la expresión [[#eq0190|(36)]] , previa resolución del movimiento del agua a través de la presa.      

En consecuencia, debe seguirse un proceso iterativo.
* 1.° Estimar una ley de variación del caudal circulante exteriormente.
* 2.° Modelar el movimiento del agua sobre la presa, determinando la carga piezométrica y el potencial a lo largo de la frontera.
* 3.° Modelar el movimiento del agua a través de la presa. Conocido el campo de velocidades puede determinarse la ley de intercambio de caudales con el exterior y la ley de variación del caudal circulante exteriormente.
* 4.° Comparar esta ley con la estimada en el primer paso.
* 5.° Si coinciden, el problema está resuelto. En caso contrario, debe modificarse la ley y reiterar el proceso.  En definitiva, vemos que la condición de frontera correspondiente a la coronación y el paramento de aguas abajo obliga, en pura teoría, a resolver el problema mediante tanteos sucesivos. Desde un punto de vista práctico, debe tenerse en cuenta que en general la altura de lámina vertiente sobre la presa será pequeña, muy inferior a las alturas de presión existentes en el dominio de filtración, por lo que la hipótesis de presión nula a lo largo del talud de aguas abajo resulta en general aceptable, salvo en el caso de que la lámina vertiente sea importante. Esta hipótesis de presión nula simplifica notablemente la modelación, ya que hace innecesaria la iteración.      

==6. Solución del problema mediante un esquema de diferencias finitas==

El método de las diferencias finitas resulta más intuitivo y conduce a una formulación más sencilla que el método de los elementos finitos. En cambio, este último tiende a ser más rápido y permite un mejor tratamiento de contornos de geometría complicada (Volker, 1969; Parkin, 1971)[[#bib0025|[5]]]  and [[#bib0020|[4]]] . Para el desarrollo de la formulación en diferentes finitas se utilizará la siguiente notación simplificada:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_x=\frac{\delta \phi }{\delta x}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_y=\frac{\delta \phi }{\delta y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}=\frac{{\delta }^2\phi }{\delta x^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{yy}=\frac{{\delta }^2\phi }{\delta y^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xy}=\frac{{\delta }^2\phi }{\delta x\delta y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)
|}

y denominaremos ''N''  = 1/''m'' , quedando la ecuación diferencial de campo de Parkin:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>({\phi }_{xx}+{\phi }_{yy})({\phi }_x^2+{\phi }_y^2)+</math><math>(N-1)({\phi }_x^2{\phi }_{xx}+2{\phi }_x{\phi }_y{\phi }_{xy}+</math><math>{\phi }_y^2{\phi }_{yy})=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)
|}

===6.1. Formulación general en un nudo interior===

Se considera una malla rectangular de lados ''b''<sub>''x''</sub>  y ''b''<sub>''y''</sub>  en las direcciones ortogonales de los ejes ''x''  e ''y''  respectivamente. Denominamos nudo interior al que está rodeado por otros 8 nudos interiores al dominio de filtración. ( [[#fig0025|fig. 5]] )

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr5.jpg|center|299px|Nudo interior.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 5.

Nudo interior.

</span>
|}

Las derivadas en un nudo interior pueden aproximarse mediante las siguientes expresiones:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_x={\left(\frac{\delta \phi }{\delta x}\right)}_O=</math><math>\frac{{\phi }_A-{\phi }_E}{2b_x}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_y={\left(\frac{\delta \phi }{\delta y}\right)}_O=</math><math>\frac{{\phi }_C-{\phi }_G}{2b_y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}={\left(\frac{{\delta }^2\phi }{\delta x^2}\right)}_O=</math><math>\frac{{\left(\frac{\delta \phi }{\delta x}\right)}_{OA}-{\left(\frac{\delta \phi }{\delta x}\right)}_{EO}}{b_x}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)
|}

Por tanto:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}=\frac{{\phi }_A-2{\phi }_O+{\phi }_E}{b_x^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 46)
|}

Por otra parte:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{yy}={\left(\frac{{\delta }^2\phi }{\delta y^2}\right)}_O=</math><math>\frac{{\left(\frac{\delta \phi }{\delta y}\right)}_{OC}-{\left(\frac{\delta \phi }{\delta y}\right)}_{GO}}{b_y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 47)
|}

y por tanto:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{yy}=\frac{{\phi }_C-2{\phi }_O+{\phi }_G}{b_y^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 48)
|}

Finalmente:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xy}={\left(\frac{{\delta }^2\phi }{\delta x\delta y}\right)}_O=</math><math>{\left[\frac{\delta }{\delta y}{\left(\frac{\delta \phi }{\delta x}\right)}_O\right]}_O=</math><math>\frac{{\left(\frac{\delta \phi }{\delta x}\right)}_C-{\left(\frac{\delta \phi }{\delta x}\right)}_G}{2\quad b_y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 49)
|}

de donde:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xy}=\frac{\frac{{\phi }_B-{\phi }_D}{2\quad b_x}-\frac{{\phi }_H-{\phi }_F}{2\quad b_x}}{2\quad b_y}=</math><math>\frac{{\phi }_B-{\phi }_D-{\phi }_H+{\phi }_F}{4\quad b_x\quad b_y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 50)
|}

Por tanto, las derivadas de la función potencial ''ϕ'' (''x'' , ''y'' ) en un nudo interior ''O'' (''i'' , ''j'' ) pueden aproximarse así:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_x(i,j)=\frac{1}{2b_x}[\phi (i,j+1)-\phi (i,j-</math><math>1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 51)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_y(i,j)=\frac{1}{2b_y}[\phi (i+1,j)-\phi (i-</math><math>1,j)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 52)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}(i,j)=\frac{1}{b_x}[\phi (i,j+1)-2\phi (i,j)+</math><math>\phi (i,j-1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 53)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{yy}(i,j)=\frac{1}{b_y^2}[\phi (i+1,j)-2\phi (i,j)+</math><math>\phi (i-1,j)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 54)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xy}(i,j)=\frac{1}{4b_xb_y}[\phi (i+1,j+1)-</math><math>\phi (i+1,j-1)-\phi (i-1,j+1)+\phi (i-1,j-1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 55)
|}

Para obtener una formulación más sencilla, llamaremos:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_x(i,j)=\frac{1}{b_x^2}[\phi (i,j+1)+\phi (i,j-1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 56)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_y(i,j)=\frac{1}{b_y^2}[\phi (i+1,j)+\phi (i-1,j)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 57)
|}

De este modo queda:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}(i,j)=K_x(i,j)-\frac{2}{b_x^2}\phi (i,j)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 58)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{yy}(i,j)=K_y(i,j)-\frac{2}{b_y^2}\phi (i,j)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 59)
|}

Por tanto, en el punto ''O'' (''i'' , ''j'' ), omitiendo los subíndices que en todos los casos son (''i'' , ''j'' ):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}+{\phi }_{yy}=K_x-\frac{2}{b_x^2}\phi +</math><math>K_y-\frac{2}{b_y^2}\phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 60)
|}

y operando:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xx}+{\phi }_{yy}=K_x+K_y-2\left(\frac{1}{b_x^2}+\right. </math><math>\left. \frac{1}{b_y^2}\right)\phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 61)
|}

Para simplificar nuevamente la formulación, llamaremos:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_{xy}=K_x+K_y</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 62)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>b_{xy}=2\left(\frac{1}{b_x^2}+\frac{1}{b_y^2}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 63)
|}

Sustituyendo en la ecuación diferencial de campo tenemos:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>(K_{xy}-b_{xy}\phi )({\phi }_x^2+{\phi }_y^2)+(N-1).\left[{\phi }_x^2\left(K_x-\right. \right. </math><math>\left. \left. \frac{2}{b_x^2}\phi \right)+2{\phi }_x{\phi }_y{\phi }_{xy}+\right. </math><math>\left. {\phi }_y^2\left(K_y-\frac{2}{b_y^2}\phi \right)\right]=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 64)
|}

Denominando:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_1={\phi }_x^2+{\phi }_y^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 65)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_2=2{\phi }_x\quad \quad {\phi }_y\quad \quad {\phi }_{xy}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 66)
|}

queda:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>(K_1\quad K_{xy}-K_1\quad b_{xy}\phi )+(N-1)\left[{\phi }_x^2\quad K_x-\right. </math><math>\left. \frac{2}{b_x^2}{\phi }_x^2\phi +K_2+{\phi }_y^2\quad K_y-\right. </math><math>\left. \frac{2}{b_y^2}{\phi }_y^2\phi \right]=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 67)
|}

Llamando ahora:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_3=K_1\quad b_{xy}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 68)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_4={\phi }_x^2\quad K_x+{\phi }_y^2\quad K_y</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 69)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_5=\frac{2{\phi }_x^2}{b_x^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 70)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_6=\frac{2{\phi }_y^2}{b_y^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 71)
|}

se obtiene:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_1K_{xy}-K_3\phi +(N-1)\quad (K_4-K_5\phi -K_6\phi )=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 72)
|}

Si llamamos:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_7=K_4+K_2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 73)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_8=K_5+K_6</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 74)
|}

tenemos:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_1\quad K_{xy}-K_3\quad \phi +(N-1)\quad (K_7-K_8\quad \phi )=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 75)
|}

de donde:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_1\quad K_{xy}-K_3\quad \phi +(N-1)\quad K_7-(N-1)\quad K_8\quad \phi =</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 76)
|}

y finalmente, denominando:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_9=K_1\quad K_{xy}+(N-1)\quad K_7</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 77)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_{10}=K_3+(N-1)\quad K_8</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 78)
|}

queda la expresión simple:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_9-K_{10}\quad \phi =0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 79)
|}

en la que puede despejarse:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{K_9}{K_{10}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 80)
|}

Deshaciendo ahora los cambios de variable realizados, salvo los de ''K''<sub>''x''</sub>  y ''K''<sub>''y''</sub> , para no complicar la expresión, se obtiene fácilmente la expresión del potencial ''ϕ'' (''i'' , ''j'' ) en un nudo interior ''O'' (''i'' , ''j'' ) en función de las derivadas ''ϕ''<sub>''x''</sub> (''i'' , ''j'' ), ''ϕ''<sub>''y''</sub> (''i'' , ''j'' ), ''ϕ''<sub>''xy''</sub> (''i'' , ''j'' ) y de los valores de ''K''<sub>''x''</sub> (''i'' , ''j'' ) y ''K''<sub>''y''</sub> (''i'' , ''j'' ), que pueden todos ponerse en función de los valores de ''ϕ''  en los nodos adyacentes, como se vio con anterioridad:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{\left({\phi }_x^2+{\phi }_y^2)(K_x+K_y)+(N-1)({\phi }_x^2K_x+{\phi }_y^2K_y+2{\phi }_x{\phi }_y{\phi }_{xy}\right)}{2\left(\frac{1}{b_x^2}+\frac{1}{b_y^2}\right)({\phi }_x^2+{\phi }_y^2)+(N-1)\left(\frac{2}{b_x^2}{\phi }_x^2+\frac{2}{b_y^2}{\phi }_y^2\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 81)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_x(i,j)=\frac{1}{b_x^2}[\phi (i,j+1)+\phi (1,j-1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 82)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_y(i,j)=\frac{1}{b_y^2}[\phi (i+1,j)+\phi (i-1,j)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 83)
|}

===6.2. Formulación en un nudo de contorno impermeable horizontal (base del espaldón)===

En un nudo situado en un contorno impermeable horizontal, como la base del espaldón, no se dispone de nudos situados bajo el mismo, pero la situación puede resolverse imponiendo la condición de que la variación del potencial de velocidad en la dirección vertical tiene que ser nula

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_y=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Esto debe cumplirse, ya que en caso contrario existiría circulación de agua en sentido vertical, lo que no es posible siendo el contorno horizontal impermeable. Imponiendo esta condición, se llega a la siguiente expresión:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{{NK}_x+K_y}{2\left(\frac{N}{b_x^2}+\frac{1}{b_y^2}\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 84)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_x(i,j)=\frac{1}{b_x^2}[\phi (i,j+1)+\phi (i,j-1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 85)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_y(i,j)=\frac{2}{b_y^2}\phi (i+1,j)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 86)
|}

===6.3. Formulación en un nudo de contorno impermeable vertical===

En este caso debe ser:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_x=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Procediendo de forma análoga obtenemos la expresión del potencial en un nudo de contorno impermeable vertical:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{K_x+N\quad K_y}{2\left(\frac{1}{b_x^2}+\frac{N}{b_y^2}\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 87)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_x(i,j)=\frac{2}{b_x^2}\phi (i,j+1)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 88)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>K_y(i,j)=\frac{1}{b_y^2}[\phi (i+1,j)+\phi (i-1,j)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 89)
|}

===6.4. Formulación en un nudo afectado por el contorno inclinado del talud===

Únicamente resulta problemática la aproximación de la derivada cruzada:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xy}={\left(\frac{{\delta }^2\phi }{\delta x\delta y}\right)}_0=</math><math>{\left[\frac{\delta }{\delta y}{\left(\frac{\delta \phi }{\delta x}\right)}_0\right]}_0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 90)
|}

De acuerdo con la notación de la [[#fig0030|figura 6]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\phi }_{xy}=\frac{{\left(\frac{\delta \phi }{\delta x}\right)}_C-{\left(\frac{\delta \phi }{\delta x}\right)}_G}{2b_y}=</math><math>\frac{{\left(\frac{\delta \phi }{\delta x}\right)}_C-\frac{{\phi }_H-{\phi }_F}{2b_x}}{2b_y}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 91)
|}

El valor de ''ϕ''<sub>''x''</sub>  en el punto ''C''  puede aproximarse de diversas formas. Se ha optado por hacer una extrapolación lineal a partir de los valores de ''ϕ''<sub>''x''</sub>  con los puntos ''D''  y ''L'' .

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>({\phi }_x)_D={\phi }_x(i+1,j-1)=\frac{1}{2b_x}[\phi (i+</math><math>1,j)-\phi (i+1,j-2)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 92)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>({\phi }_x)_L=\frac{1}{2\left(\frac{b_x}{2}\right)}[\phi (i+</math><math>1,j)-\phi (i+1,j-1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 93)
|}

Extrapolando linealmente:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>({\phi }_x)_C=({\phi }_x)_D+2[({\phi }_x)_L-({\phi }_x)_D]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 94)
|}

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr6.jpg|center|336px|Nudo afectado por el contorno inclinado del talud.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 6.

Nudo afectado por el contorno inclinado del talud.

</span>
|}

===6.5. Punto singular intersección de dos contornos impermeables, uno horizontal y el otro vertical===

En este punto singular el potencial de velocidad se puede aproximar considerando el valor medio de las extrapolaciones lineales realizadas a partir de los dos nudos más próximos en fila y columna, es decir, en el contorno impermeable horizontal y vertical cuya intersección origina el punto singular. Considerando dicho punto singular como origen de filas y columnas, tendrá ''i''  = 0 y ''j''  = 0. La extrapolación a partir de los dos nodos más próximos de la fila 1, contorno horizontal impermeable, será:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi (1,1)^\ast =\phi (1,3)+2[(\phi (1,2)-\phi (1,3)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 95)
|}

y la extrapolación a partir de los dos nodos más próximos de la columna 1, contorno vertical impermeable:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi (1,1)^{\ast \ast }=\phi (3,1)+2[\phi (2,1)-\phi (3,1)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 96)
|}

Por tanto:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi (1,1)=\frac{1}{2}[\phi (1,1)^\ast +\phi (1,1)^{\ast \ast }]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 97)
|}

Esta formulación en diferencias finitas se ha implementado en el código que se ha denominado Movimiento No Lineal en Espaldón de Escollera. Red de filtración (MNLEE.RDF).      

==7. Resultados del modelo numérico==

Una vez resuelto el problema definido por la ecuación diferencial de campo y las condiciones de contorno (mediante el código MNLEE.RDF), el potencial de velocidad (''ϕ'' ) es conocido en todos los nudos del dominio de filtración. A partir de ellos puede obtenerse también en cada nudo la velocidad de filtración y el gradiente hidráulico, que definen sendos campos vectoriales, y el campo escalar de presiones del agua. También puede obtenerse el caudal filtrado a través de la escollera.      

===7.1. Presión del agua===

El potencial de velocidad en un punto del dominio de filtración es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{h}{c}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 98)
|}

siendo «''h'' »la carga hidráulica o carga piezométrica, ya que la carga de velocidad es despreciable frente a la suma de altura y presión.

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>h=z+\frac{p}{{\gamma }_w}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 99)
|}

y ''c''  la constante de la fórmula de resistencia:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i=c\quad V^m</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 100)
|}

Por tanto, conocido ''ϕ'' (''i'' , ''j'' ) en cada nudo del dominio, una vez resuelto el problema, de forma inmediata puede obtenerse la carga ''h'' (''i'' , ''j'' ) correspondiente.

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>h(i,j)=c\quad \phi (i,j)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 101)
|}

y la presión del agua

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>p(i,j)={\gamma }_w[h(i,j)-z(i,j)]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 102)
|}

Si bien la presión es fluctuante en las proximidades del talud, debido a la turbulencia, la fluctuación se atenúa rápidamente hacia el interior del medio poroso (Detert et al., 2010)[[#bib0030|[6]]] .      

===7.2. Velocidad de filtración===

Como se vio al presentar la ecuación de campo de Parkin, las componentes del vector velocidad de filtración son:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u=\frac{\delta \phi }{\delta x}{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 103)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>v=\frac{\delta \phi }{\delta y}{\left[{\left(\frac{\delta \phi }{\delta x}\right)}^2+{\left(\frac{\delta \phi }{\delta y}\right)}^2\right]}^{\frac{1}{2}\left(\frac{1}{m}-1\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 104)
|}

A partir de estas ecuaciones puede obtenerse el vector velocidad de filtración en cada nudo del dominio. Su módulo será:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>V=\sqrt{u^2+v^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 105)
|}

y el ángulo de avance respecto a la horizontal

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\alpha =arctan\left(\frac{v}{u}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 106)
|}

con sentido antihorario positivo.      

===7.3. Gradiente hidráulico===

El vector gradiente hidráulico en un punto lleva la dirección de la velocidad de filtración, o sea, es tangente a la línea de filtración que pasa por ese punto, y su módulo puede obtenerse directamente a partir de su propia definición o bien aplicando la fórmula de resistencia, una vez conocida la velocidad de filtración en el punto considerado. Por la propia definición de gradiente hidráulico:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{i}=\boldsymbol{grad}\mbox{ }h</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 107)
|}

Teniendo en cuenta que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\phi =\frac{h}{c}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 108)
|}

donde ''c''  es una constante, tenemos:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{i}=c\mbox{ }\boldsymbol{grad}\mbox{ }\phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 109)
|}

Las coordenadas de <math display="inline">\boldsymbol{grad}\mbox{ }\phi </math> son:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{\delta \phi }{\delta x}\quad \quad ,\quad \quad \frac{\delta \phi }{\delta y}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 110)
|}

Por tanto, el módulo del gradiente hidráulico puede obtenerse en el nudo (''i'' , ''j'' ):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i(i,j)=c\quad \quad \sqrt{{\left[\frac{\delta \phi }{\delta x}(i,j)\right]}^2+{\left[\frac{\delta \phi }{\delta y}(i,j)\right]}^2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 111)
|}

En cada nudo los valores de las derivadas pueden obtenerse mediante las expresiones deducidas con anterioridad. También puede aplicarse la fórmula de resistencia, una vez calculada la velocidad en cada nudo como se vio anteriormente. El gradiente hidráulico tendrá el siguiente módulo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i(i,j)=c\quad \quad [V(i,j)]^m</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 112)
|}

Obviamente, el valor obtenido por medio de cualquiera de las dos expresiones debe ser idéntico en todos los nudos, como se ha comprobado mediante el código MNLEE.RDF.

===7.4. Caudales filtrados===

Admitiendo que el espaldón de aguas abajo está totalmente saturado por el agua que vierte sobre la coronación de la presa, el caudal que circula a través de distintas secciones verticales es distinto, ya que entre cada 2 secciones se produce una aportación o detracción de agua a través del talud ([[#fig0035|fig. 7]] ).

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr7.jpg|center|339px|Intercambio de caudales entre el interior y el exterior del espaldón de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 7.

Intercambio de caudales entre el interior y el exterior del espaldón de escollera.

</span>
|}

Denominaremos: ''Q''<sub>''c''</sub> (''j'' ) el caudal que circula a través de la sección definida por la columna ''j'' , ''Q''<sub>''f''</sub> (''j'' ) caudal que se filtra a través del talud entre las secciones definidas por las columnas «''j'' »y «''j''  + 1». El caudal ''Q''<sub>''c''</sub> (''j'' ), y análogamente ''Q''<sub>''c''</sub> (''j''  + 1), puede obtenerse por suma de los caudales circulantes entre cada dos nudos. Conocida la componente horizontal de la velocidad en cada nudo, puede considerarse una velocidad media en el interespacio entre cada dos nudos consecutivos en la misma vertical. Por tanto:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q_c(j)=\sum_{i=1}^{i=n_j}\frac{1}{2}[u(i,j)+u(i+1,j)]\quad \quad b_y</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 113)
|}

siendo ''n''<sub>''j''</sub>  el número de nudos de la columna ''j'' .      

Por diferencia entre los caudales circulantes entre cada dos secciones verticales, definidas por 2 columnas consecutivas, puede obtenerse el caudal que se filtra a través del talud, con sentido positivo, de entrada, o negativo, de salida, entre las dos secciones consideradas:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q_f(j)=Q_c(j+1)-Q_c(j)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 114)
|}

Es claro que la suma de todos los caudales ''Q''<sub>''f''</sub> (''j'' ) debe ser nula, ya que el caudal que entra al cuerpo de presa por la parte superior sale por la parte inferior. Esto se ha comprobado al utilizar el código MNLEE.RDF.      

Por otra parte, la suma de todos los caudales ''Q''<sub>''f''</sub> (''j'' ) positivos define el caudal máximo que circula por el interior del espaldón, que se denominará caudal de saturación y se denotará por ''Q''<sub>''f''</sub> . Por tanto, siguiendo el proceso expuesto puede definirse la ley de caudales entrantes y salientes a través del talud sometido a vertido, que se denominará ley de intercambio de caudales así como el caudal que circula por el interior del espaldón, que es variable al considerar sucesivas secciones debido al intercambio de caudal que se produce a través del talud, y particularmente se podrá calcular el valor del caudal de saturación. También es posible, observando la sección en que el caudal ''Q''<sub>''f''</sub> (''j'' ) pasa de ser positivo a ser negativo, conocer la longitud de la zona de infiltración del talud, a través de la cual el agua entra en el espaldón y la longitud de la zona de emergencia, a través de la cual el agua sale del espaldón.      

==8. Código MNLEE.RDS==

Como se indicó anteriormente, la formulación definida permitió implementar un código para realizar el cálculo de la red de filtración y las presiones asociadas y caudales circulantes considerando la no linealidad del movimiento. Las principales hipótesis consideradas en dicho programa se resumen en:
* Geometría definida en la [[#fig0040|figura 8]]

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr8.jpg|center|323px|Definición geométrica de la sección del espaldón en el código MNLEE.RDS.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 8.

Definición geométrica de la sección del espaldón en el código MNLEE.RDS.

</span>
|}
* Espaldón completamente saturado.
* Situación estacionaria.
* Movimiento no lineal (fórmula de resistencia exponencial y ecuación de campo de Parkin).
* Lámina de agua sobre el talud despreciable.

La determinación del potencial de velocidad se realiza mediante un procedimiento iterativo. Se ha definido inicialmente un valor máximo del error admisible en la determinación del potencial de velocidad en cada punto de la malla. Después de realizar el cálculo del potencial en cada nudo de dicha malla, se ha comprobado el error cometido. En el caso de que el error sea mayor que el máximo admisible se realiza una corrección del valor del potencial para volver a realizar el cálculo. El problema se considera resuelto cuando el potencial de velocidad está ajustado en todos los puntos de la malla.

En el caso de que la lámina de agua no fuera despreciable con relación al campo de potenciales existentes en el interior del dominio de filtración, el diagrama de flujo habría de completarse introduciendo un proceso iterativo para utilizar la ley de intercambio de caudales como elemento de ajuste.

==9. Campos de presiones, velocidades y gradientes hidráulicos en una presa de escollera sometida a sobrevertido==

El código MNLEE.RDF de cálculo numérico desarrollado ha permitido modelar el movimiento del agua a través del espaldón de aguas abajo de diversos casos de presas de escollera obteniendo como resultado los campos de presiones, velocidades y gradientes hidráulicos, así como cuantificar el caudal de sobrevertido necesario para saturar completamente el referido espaldón.

En la [[#fig0045|figura 9]]  se representa el campo de presiones intersticiales en el espaldón de una presa de 50 ''m  ''  de altura y talud 1,5:1 sometido a la circulación de agua como consecuencia de un vertido sobre la coronación de la presa. La presión intersticial se expresa en metros de columna de agua, tanto en esta figura como en las demás. La modelación numérica del fenómeno de filtración se ha realizado considerando una relación no lineal entre el gradiente hidráulico y la velocidad de filtración dada por la ley de resistencia <math display="inline">i=7,13\quad \quad {10}^{-2}\quad \quad V_d^{1.85}</math> . Se observa que las líneas isobaras tienden, en su recorrido de aguas arriba hacia aguas abajo, a colocarse progresivamente cada vez más paralelas al talud de la presa. El paralelismo ocurre cuando la presión del agua es la hidrostática, valor máximo que puede alcanzar la presión del agua en un movimiento de lámina libre. Por tanto, el patrón de isobaras observado indica que las presiones se aproximan a ese máximo posible para cada punto, en función de su situación, en la zona del pie de presa. Esto deberá tenerse en cuenta a la hora de analizar la estabilidad de la presa frente al deslizamiento en masa en situación de sobrevertido.

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr9.jpg|center|357px|Campo de presiones intersticiales representado mediante líneas isobaras.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 9.

Campo de presiones intersticiales representado mediante líneas isobaras.

</span>
|}

El patrón de isolíneas del módulo del gradiente hidráulico indica que los valores máximos se producen en las proximidades del talud, especialmente en las zonas de coronación y pie de presa, mientras que los gradientes hidráulicos mínimos se dan en la zona más interior del espaldón ([[#fig0050|fig. 10]] ). El gradiente hidráulico es un factor crítico en el fenómeno de arrastre de partículas. Hay que destacar que en la zona de coronación el agua entra en la presa, y por tanto la fuerza de arrastre que depende del gradiente se dirige hacia el interior del espaldón, por lo que esta zona no resulta problemática en lo que se refiere al arrastre de partículas. En cambio, en la zona de pie de presa la dirección del agua es hacia el exterior, y por tanto un elevado gradiente hidráulico contribuye, a través de una mayor fuerza de arrastre, a aumentar el riesgo de pérdida de partículas.

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr10.jpg|center|357px|Campo de gradientes hidráulicos expresado mediante isolíneas.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 10.

Campo de gradientes hidráulicos expresado mediante isolíneas.

</span>
|}

Las isolíneas del módulo de la velocidad de filtración presentan un patrón similar al de las isolíneas del módulo del gradiente hidráulico, como cabía esperar ([[#fig0055|fig. 11]] ). En la [[#fig0060|figura 12]]  puede observarse el módulo y la dirección de la velocidad de filtración en cada nodo del espaldón de una presa de 50 ''m  ''  de altura y talud 2:1, con ley de resistencia de la escollera <math display="inline">i=0,40\quad \quad V_d^{1.85}</math> . De ello se induce un patrón de líneas de corriente semejante al obtenido mediante modelación física en laboratorio (Gerodetti, 1981)[[#bib0035|[7]]] , con entrada de agua por la parte superior del talud y salida por la zona inferior.

<span id='fig0055'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr11.jpg|center|357px|Campo del módulo de la velocidad de filtración, expresado mediante isolíneas.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 11.

Campo del módulo de la velocidad de filtración, expresado mediante isolíneas.

</span>
|}

<span id='fig0060'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr12.jpg|center|358px|Módulo y dirección de la velocidad de filtración.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 12.

Módulo y dirección de la velocidad de filtración.

</span>
|}

Puede definirse el caudal de saturación como el caudal necesario para saturar completamente la escollera del espaldón. Este caudal depende principalmente de la permeabilidad del material que constituye el espaldón y de la altura de la presa, siendo la influencia del talud mínima, despreciable a efectos prácticos, según ha podido establecerse mediante la realización de un estudio paramétrico. Para una ley de resistencia <math display="inline">i=7,13\quad \quad {10}^{-2}\quad \quad V_d^{1.85}</math>  el caudal de saturación unitario, por metro de anchura de presa, oscila entre algo menos de 0, 5 ''m''<sup>2</sup> /''s''  para una presa de 25 ''m''  de altura y algo más de 2 ''m''<sup>2</sup> /''s''  para una presa con una altura de 150 ''m''  ( [[#fig0065|fig. 13]] ). Hay que tener en cuenta que estos valores son unitarios, de modo que los caudales de saturación totales pueden llegar a ser importantes en presas con elevadas longitudes de coronación y anchura de valle. El cálculo realizado es bidimensional, y por tanto no se ha tenido en cuenta la convergencia de las laderas, por lo que los resultados son aplicables esencialmente a presas de gran longitud, en las que el efecto de la convergencia de laderas es de segundo orden.

<span id='fig0065'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr13.jpg|center|343px|Caudal de saturación del espaldón de aguas abajo, según la altura de presa y el ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 13.

Caudal de saturación del espaldón de aguas abajo, según la altura de presa y el talud.

</span>
|}

==10. Efecto de la no linealidad del movimiento del agua a través de la escollera==

La mayoría de los códigos comerciales para cálculo de redes de filtración implementan la ley de Darcy como fórmula de resistencia al avance del agua, es decir, admiten la linealidad del movimiento. Dado que en materiales granulares gruesos como los que constituyen generalmente los espaldones de una presa de escollera la ley de resistencia no es lineal, cabe preguntarse cuál es la influencia práctica de la no linealidad.

Se han modelado mediante el código MNLEE.RDS dos casos de presas de 50 ''m''  de altura con taludes de 1,5:1 y 3:1, considerando una relación entre el gradiente hidráulico y la velocidad no lineal de tipo exponencial con exponente 1,85. Los mismos casos se han modelado mediante un código comercial (SEEP/W), admitiendo la ley de Darcy. A partir de los resultados de presiones y gradientes hidráulicos obtenidos, se han determinado las diferencias, que vienen a representar el error de linealización, es decir, el error cometido como consecuencia de asumir una ley de Darcy cuando en realidad la ley de resistencia no es lineal. Los resultados pueden verse en las figuras  [[#fig0070|14]]  y [[#fig0075|15]] .

<span id='fig0070'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr14.jpg|center|358px|Error de linealización, en presiones.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 14.

Error de linealización, en presiones.

</span>
|}

<span id='fig0075'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_163438791-1-s2.0-S0213131512000144-gr15.jpg|center|358px|Error de linealización, en gradientes hidráulicos.]]


|-
| <span style="text-align: center; font-size: 75%;">

Fig. 15.

Error de linealización, en gradientes hidráulicos.

</span>
|}

Respecto de las presiones intersticiales, se observa que una línea aproximadamente ortogonal al talud separa 2 zonas en las que el error de linealización es de signo contrario: las presiones obtenidas mediante el modelo lineal son menores en la zona de pie de presa y mayores en la zona de coronación que las obtenidas mediante el modelo no lineal. Los errores máximos se producen en zonas internas del espaldón, alejadas del talud, donde los errores son mínimos. Dado que la estabilidad de la presa depende esencialmente de las presiones intersticiales en las proximidades del talud, con los resultados obtenidos puede decirse que las presiones obtenidas mediante un modelo lineal tienen una precisión suficiente a efectos prácticos para estudiar a partir de ellas la estabilidad de la presa. Por otra parte, se aprecia que con el talud más suave el error de linealización es menor.

Respecto de los gradientes hidráulicos, se aprecia que el error de linealización es por defecto en las zonas de coronación y pie de presa, mientras que en una amplia zona central el error es por exceso. En la zona crítica de pie de presa el error de linealización no supera el valor 0,02 en ninguno de los casos, lo que supone un error moderado respecto de los gradientes hidráulicos en pie de presa de 0,67 para el caso de talud 1,5:1 y de 0,33 para el caso de talud 3:1.

==11. Conclusiones==

El movimiento del agua sobre y a través del espaldón de aguas abajo de una presa de escollera sometida a sobrevertido tiene algunas características particulares: a) el movimiento del agua a través de la escollera es no lineal en el sentido de que la relación entre la velocidad de filtración y el gradiente hidráulico no es lineal; b) el movimiento del agua externo, sobre el talud de escollera, es espacialmente variado, con caudales decrecientes en la zona superior, de infiltración, y crecientes en la zona inferior, de emergencia, y c) los movimientos de filtración a través de la escollera y exterior, sobre el talud de aguas abajo, están acoplados. La ley de intercambio de caudales, que define los caudales unitarios que se infiltran hacia el interior del espaldón o que emergen al exterior por el mismo, puede utilizarse como elemento de ajuste para resolver el problema acoplado de forma iterativa. Si los calados de agua sobre el talud son suficientemente reducidos frente al campo de presiones interior, el problema de filtración puede resolverse de forma independiente obteniendo como resultado los campos de presiones, velocidades y gradientes en todo el espaldón de escollera, y también la ley de intercambio de caudales. Esto puede hacerse mediante el esquema de diferencias finitas que aquí se presenta, y que se ha implementado en el código MNLEE.RDS. La resolución del problema de filtración y, por tanto, el conocimiento de los campos de presiones, velocidades de filtración y gradientes hidráulicos, proporciona los datos necesarios para abordar el cálculo de la estabilidad de la presa y del problema de arrastre de partículas en una presa de escollera sometida a sobrevertido, abriendo la posibilidad de estudiar la seguridad de este tipo de presas en tal situación. También es posible acometer el desarrollo de criterios de diseño para que las presas de escollera sean seguras hasta un determinado caudal de sobrevertido.

==Agradecimientos==

Al Ministerio de Ciencia e Innovación, como entidad financiadora del proyecto de investigación del Plan Nacional de I+D 2007-2011 denominado «Caracterización de la rotura de las presas de escollera por sobrevertido y desarrollo de criterios para evaluar la seguridad del conjunto presa-área afectada durante una avenida», (XPRES) con código de identificación BIA2007-68120-C03-02.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] G.H. George, D. Hansen; Conversion between quadratic and power law for non-Darcy flow, ASCE; J. Hydraul. Eng., 118 (5) (1992), pp. 792–797</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] A. Parkin, D.H. Trollope, J.D. Lawson; Rockfill structures subject to water flow, ASCE; J. Soil. Mech. Foundations. Div. Asce., 92 (6) (1966), pp. 135–151</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] F. Engelund; On laminar and turbulent flow of groundwater in homogeneous sand; Trans. Dan. Acad. Tech. Sci., 3 (1953), pp. 1–105</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] A. Parkin; Field solutions for turbulent seepage flow, ASCE; J. Soil Mech. Found. Div., 97 (1) (1971), pp. 209–218</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] R.E. Volker; Nonlinear flow in porous media by finite elements, ASCE; J. Hydraulics Div., 95 (6) (1969), pp. 2093–2114</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] M. Detert,  ''et al.''; Laboratory measurements on turbulent pressure fluctuations in and above gravel beds, ASCE; J. Hydraulic Eng., 136 (10) (2010), pp. 779–789</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] M. Gerodetti; Model studies of an overtopped rockfill dam; Water Power & Dam Construct., 33 (9) (1981), pp. 25–31</li>
</ol>