En otros tutoriales hemos aprendido a utilizar ArcGIS para realizar algunos análisis hidrológicos como la obtención de Parámetros de cuenca con ArcGIS o la generación de una Red hidrológica de órdenes. Además, tal y como vimos en el tutorial de Análisis de erosión con ArcGIS, uno de los factores que intervienen en la Ecuación Universal de la pérdida de suelo es la erosividad de la lluvia. En esta ocasión vamos a calcular el índice de erosividad de la lluvia con ArcGIS.
Cálculo del índice de erosividad de la lluvia con ArcGIS
La Erosividad de la lluvia se puede definir como su agresividad sobre el suelo. Representa la energía con que las gotas de lluvia que caen a una determinada intensidad pueden romper los agregados superficiales en partículas de tamaño transportable.
La erosividad es función de la energía cinética (E) y de la intensidad de la precipitación (I). La primera, relaciona el tamaño y la velocidad de caída de las gotas y la segunda relaciona la cantidad con la duración. Por tanto, se trata de un parámetro complejo de calcular pues relaciona cuatro propiedades (cantidad y duración, diámetro y velocidad de las gotas de lluvia).
El índice propuesto por la ecuación universal de la pérdida de suelo (USLE) para cuantificar este efecto erosivo de las lluvias debido a su intensidad presenta la siguiente expresión:
donde R representa el Índice de Erosividad de la lluvia, n es el número de aguaceros del periodo considerado, I30 es la intensidad máxima en 30 minutos para cada evento [mm/h] y E es la energía cinética del aguacero que se estima mediante las siguientes fórmulas [MJ/ha].
- Si I <= 76 mm/h –> E = 0.119 + 0.0873 logI
- Si I > 76 mm/h –> E = 0.283
Para realizar estos cálculos se debe disponer de una información muy detallada de las precipitaciones ya que se debe contar con un registro continuo de las variaciones de intensidad de la lluvia durante los diferentes aguaceros. Por este motivo, en un intento de facilitar dicho cálculo, el antiguo Instituto para la Conservación de la Naturaleza (I.C.O.N.A.) a través de su publicación «Agresividad de la lluvia en España» (MAPA , 1988) calculó para la Península Ibérica unas ecuaciones para sustituir a la R = EI30 de WISCHMEIER Y SMITH, de manera que con unos datos más sencillos, se obtienen resultados similares. Se trata de tres ecuaciones diferentes que se utilizan en función de la zona geográfica donde se realiza el estudio.
Pero además del I.C.O.N.A. diversos autores han intentado relacionar el factor R con parámetros más fáciles de obtener y calcular. Con este propósito, FOURNIER (1960) estableció el índice de agresividad climática o Índice de Fournier (IF) cuyo cálculo se realiza a partir de los datos pluviométricos de estaciones meteorológicas representativas, según la siguiente ecuación:
donde IF es el Índice de Fournier, Pmáx es la precipitación media correspondiente al mes más lluvioso [mm] y P es la precipitación media anual [mm].
Sin embargo, éste índice únicamente considera el mes de mayor precipitación por lo que debe ser utilizado con cautela, especialmente en aquellas zonas que presentan un régimen pluvial con más de un pico mensual de precipitación, o donde los valores pluviométricos son en general elevados.
Precisamente por este motivo ARNOLDUS (1978) propuso una corrección del Índice de Fournier para considerar no sólo la precipitación mensual del mes más húmedo, sino también la del resto de los meses. Este Índice de Modificado de Fournier (IMF) caracteriza la agresividad de la precipitación y se calcula de la siguiente forma:
donde pi es la precipitación de cada mes [mm] y Pt es la precipitación media anual [mm].
En cualquier caso, la gran mayoría de las veces es muy necesario realizar ajustes de regresión entre el factor R calculado por el método USLE y el correspondiente IF o IMF. Así, en la literatura es posible encontrar diferentes estudios en los que se aplican distintos métodos que correlacionan R con el IF o el IMF mediante ecuaciones diferentes en función de la región donde se esté realizando el estudio.
Uno de ellos es el realizado por Arnoldus (1977) que estableció un mapa de índices de erosividad R para Marruecos empleando el IMF en la forma:
R = 0.264 * IMF ^ 1.50
O RODRIGUEZ ET AL. (2004) que emplearon la ecuación genérica del I.C.O.N.A. para estimar R en la zona de Badajoz:
CÁLCULO DEL ÍNDICE DE EROSIVIDAD DE LA LLUVIA CON ARCGIS
Para calcular el índice de erosividad de la lluvia (factor R), partimos de un Excel con los datos de precipitaciones medias mensuales para 279 estaciones. El Excel recoge el código de cada estación y las precipitaciones desde enero hasta diciembre.
Por otro lado, tenemos un shape con la ubicación de cada estación, cuyo código es el mismo que en el Excel así que lo primero que haremos será un Join para obtener un shape con todos los datos: capa estaciones (botón derecho) < Joins & Relates < Join.
La unión se hace siempre mediante un campo común entre la capa y la tabla (archivo) que se van a unir. En nuestro caso, la capa contiene el código de cada estación, y el archivo ‘Precip_mens.xls’ contiene los valores de precipitación, por lo tanto, lógicamente el campo común será ‘Código’.
Procedemos a indicar cada uno de los parámetros que nos pide el programa: primero elegimos el campo de la capa en el que se basa la unión (seleccionamos ‘Código’ con el desplegable). En segundo lugar indicamos la tabla que queremos unir; si ya la hemos cargado en el proyecto podemos seleccionarla con el desplegable; en caso contrario buscamos el archivo correspondiente con el icono de la carpeta; una vez realizado tendremos todos los datos en nuestra capa de estaciones:
Para poder calcular la erosividad de la lluvia, debemos disponer de las precipitaciones en formato raster para cada mes. Para obtener una información de este tipo debemos realizar una interpolación.
Existen diferentes métodos de interpolación y nosotros vamos a usar el método de la Ponderación de Distancia Inversa (IDW).
Mediante este método se calcula el valor de cada celda como la media ponderada de los valores del entorno en función del inverso de la distancia, por lo que se asume que los puntos más cercanos tendrán más influencia. El método asigna pesos a los datos del entorno en función inversa de la distancia que los separa de la celda que se estima. Esto implica que cuanto más cerca está un punto del centro de la celda, más influencia o peso tiene en el valor que adquiere.
Para determinar el grado de influencia que tendrán los puntos en función de la distancia, ésta se eleva a una potencia determinada (el peso de cada punto es proporcional al inverso de la distancia elevado a esa potencia).
El valor de la potencia es configurable: si queremos otorgar mayor importancia a los puntos más cercanos, asignaremos valores de potencia altos y obtendremos como resultado una superficie más detallada; por el contrario, al asignar valores de potencia más bajos indicamos mayor influencia de los puntos más alejados dando como resultado una superficie más suavizada (por lo general se suele utilizar un valor de potencia de 2).
Otro de los parámetros que podemos configurar es el Radio de búsqueda. Mediante este parámetro es posible limitar los puntos que consideraremos para calcular el valor de cada celda de salida. Este parámetro puede resultar de gran ayuda puesto que en muchas ocasiones los puntos más alejados de la celda que se está calculando tendrán escasa influencia en su valor por lo que sería innecesario incluirlos en el cálculo.
ArcGIS presenta dos opciones diferentes para establecer el radio de búsqueda y en función del que elijamos tendremos que indicar la cantidad de puntos a considerar (un número concreto o un mínimo) y/o una distancia (específica o máxima).
Opcionalmente, se pueden introducir líneas de corte para limitar la búsqueda, es lo que se conoce como «barreras».
Vamos a realizar el cálculo correspondiente seleccionando este método en la caja de herramientas de análisis espacial de interpolación ArcToolbox < Spatial Analyst Tools < Interpolation < IDW.
Seleccionamos la capa que contiene las estaciones como puntos de entrada y como ‘Z_value_field’ debemos indicar el campo de la tabla de atributos que queremos interpolar. Comenzamos interpolando los datos del mes de enero. Para este estudio nosotros vamos a dejar los valores potencia y radio de búsqueda que aparecen por defecto (Potencia ‘Power’ = 2, radio de búsqueda ‘Search radius’ = Variable y puntos a considerar ‘Number of Points’ = 12). Además vamos a trabajar con un tamaño de celda de 100 metros (‘Output cell size’ = 100).
Tras ejecutar la herramienta hemos obtenido un raster de 100 metros de tamaño de celda que representa el valor de precipitación de enero en cada ubicación.
De esta manera, realizamos la misma operación con los datos de precipitación de cada mes.
Una vez realizado contamos 12 raster que contienen la precipitación media mensual de manera que cada uno de los raster se corresponde con un mes, desde enero a diciembre.
Una vez que tenemos cargados los 12 raster en el proyecto comenzamos a calcular el índice de erosividad de la lluvia que, tal y como hemos explicado, se basará en la fórmula R = 2,56* IMF^1,065
Esto implica que tendremos que calcular el Índice Modificado de Fournier:
donde pi es la precipitación de cada mes [mm] y Pt es la precipitación media anual [mm].
INDICE MODIFICADO DE FORUNIER
Cada uno de los raster se corresponde por tanto con cada uno de los 12 valores de pi por lo que a partir de ellos podemos obtener el raster que representa el valor de la precipitación media anual (Pt). Este valor será la suma de todos los valores mensuales de precipitación.
Para calcularlo utilizamos la Calculadora raster que es una herramienta que permite realizar distintas operaciones matemáticas entre capas raster ArcToolbox < Spatial Analyst Tools < Map Algebra < Raster Calculator.
La Calculadora raster nos muestra todas y cada una de las capas disponibles en el proyecto actual con las que podemos realizar las operaciones.
Lo que debemos hacer es sumar los 12 raster de precipitación mensual; para ello seleccionamos uno haciendo doble click sobre la capa, pulsamos sobre el símbolo «+», seleccionamos la segunda capa a sumar y así sucesivamente hasta obtener la expresión que se muestra en la siguiente imagen. Por último, tendremos que asignar un nombre a la nueva capa que se va a generar
Obtenemos un raster que contiene los valores de precipitación media anual, que en este caso, van desde un mínimo de 328 mm hasta los 1762 mm.
Por tanto ya tenemos los valores de pi y de Pt así que ya podemos calcular IMF según la fórmula correspondiente.
Utilizamos de nuevo la Calculadora raster; en este caso tenemos que elevar al cuadrado cada uno de los valores de pi y para ello utilizamos la opción denominada «Square» que podemos seleccionar en la parte derecha de la ventana. Así iremos escribiendo la expresión correspondiente y obtendremos un IMF cuyo rango de valores va desde 35 a 213 mm.
INDICE DE EROSIVIDAD DE LA LLUVIA
Por último calculamos el valor de R según la expresión que hemos explicado anteriormente; abrimos la Calculadora raster una vez más.
Para aplicar esta fórmula tendremos que utilizar la expresión «Power» que permite elevar cualquier número a una potencia; al seleccionarla nos aparecerá un paréntesis donde el primer número será el IMF y el segundo será valor de la potencia a la que lo queremos elevar.
Finalmente obtenemos el raster que contiene los valores del índice de erosividad de la lluvia (R) que se mueve en el rango de 113 y 773.
El mayor valor de erosividad de la lluvia que se alcanza es de unos 773 [MJ.mm/ha.h] presentando un mínimo de aproximadamente 113 [MJ.mm/ha.h].
Observamos que los valores más elevados se concentran en la zona occidental de la región y puntualmente en una zona concreta hacia el este lo que contrasta con los valores más bajos que se expanden en el eje norte-sur en la franja intermedia de la región.
En base a este factor, cabe esperar que el mayor nivel de erosión se produzca en las zonas con mayor índice de erosividad de la lluvia, es decir, la zona oeste.
