UNIVERSIDADE ESTADUAL PAULISTA - UNESP CÂMPUS DE JABOTICABAL MODELAGEM DA PERDA DE FÓSFORO POR EROSÃO HÍDRICA Vera Lucia da Silva Farias Bióloga 2016 UNIVERSIDADE ESTADUAL PAULISTA – UNESP CÂMPUS DE JABOTICABAL MODELAGEM DA PERDA DE FÓSFORO POR EROSÃO HÍDRICA Autora: Vera Lucia da Silva Farias Orientador: Prof. Dr. Marcilio Vieira Martins Filho Coorientador: Prof. Dr. Millôr Godoy Sabará Tese apresentada à Faculdade de Ciências Agrárias e Veterinárias - UNESP, Câmpus de Jaboticabal, como parte das exigências para a obtenção do título de Doutora em Agronomia (Ciência do Solo) 2016 Farias, Vera Lucia da Silva F224m Modelagem da perda de fósforo por erosão hídrica / Vera Lucia da Silva Farias. – – Jaboticabal, 2016 v, 61 p. : il. ; 29 cm Tese (doutorado) - Universidade Estadual Paulista, Faculdade de Ciências Agrárias e Veterinárias, 2016 Orientador: Marcilio Vieira Martins Filho Banca examinadora: Daniela Tolêdo de Paula, Diego Silva Siqueira, Ludmila de Freitas, Gener Tadeu Pereira Bibliografia 1. Taxa de enriquecimento. 2. Modelo WEPP. 3. Variabilidade espacial. I. Título. II. Jaboticabal-Faculdade de Ciências Agrárias e Veterinárias. CDU 631.459 Ficha catalográfica elaborada pela Seção Técnica de Aquisição e Tratamento da informação – Serviço Técnico de Biblioteca e Documentação - UNESP, Câmpus de Jaboticabal. E-mail: verlucbio@yahoo.com.br mailto:verlucbio@yahoo.com.br DADOS CURRICULARES DO AUTOR VERA LÚCIA DA SILVA FARIAS - Filha de José Pedro da Silva e Lúcia Nonato da Silva. Nasceu no município de Bom Despacho - MG, em 08-06-1970. Cursou Licenciatura em Biologia, no Centro Universitário do Triangulo Mineiro (UNITRI), e concluiu em 2005. Em 2008, concluiu o Latu Sensu em Análises Clínicas - Avanços Diagnósticos, na área de concentração de Ciências Biológicas, da Faculdade de Medicina de São José do Rio Preto (FAMERP). Ingressou no mestrado em Agronomia, na área de Ciência do Solo, em março de 2011, na Universidade Estadual Paulista (UNESP/FCAV), em Jaboticabal- SP. Iniciou o curso de Doutorado em Agronomia, na área de Ciência do Solo, em março de 2013, na Universidade Estadual Paulista (UNESP/FCAV), em Jaboticabal- SP. É docente na Faculdade Frutal (FAF), desde 2008 até a presente data e na Universidade do Estado de Minas Gerais (UEMG) de 2009 até o momento. Professora efetiva de Ciências/Biologia no Estado de Minas Gerais. Se Eu fechar os céus, e não houver chuva; ou se ordenar aos gafanhotos que consumam a terra; ou se enviar a peste entre o meu povo; E se o meu povo, que se chama pelo Meu nome, se humilhar, e orar, e buscar a minha face e se converter dos seus maus caminhos, então Eu ouvirei dos céus, e perdoarei os seus pecados, e sararei a sua terra. II Crônicas 7:13-14 Ao maior tesouro da minha vida: minha família. DEDICO A todos aqueles que se sentem felizes por mais uma vitória em minha vida. OFEREÇO AGRADECIMENTOS Agradeço, primeiramente, a Deus, por me abençoar e capacitar diante de todos os passos de minha vida. À minha família, meu esposo Edvaldo Farias de Souza, meus filhos: Pedro e Luiza Farias da Silva Souza, por estarem sempre presentes na minha vida, me incentivando, dando forças e transmitindo todo o amor deste mundo. Amo vocês. Á Primeira Igreja Batista em Frutal, por estar presente e em oração na minha vida. À Universidade Estadual Paulista, UNESP, pelas grandes oportunidades de diálogo entre diversas áreas, conhecimentos, pessoas, que contribuíram para minha formação, não apenas acadêmica, mas como experiência de vida e formação como “ser” humano em todas suas esferas. À Faculdade Frutal pela oportunidade de desenvolvimento profissional. À Universidade do Estado de Minas Gerais, UEMG, também, pela oportunidade de desenvolvimento profissional, Ao Professor Marcílio Vieira Martins Filho, por sua orientação, e se disponibilizar sempre, tão atencioso e dedicado a responder minhas inquietudes e dúvidas que me afligiram durante um bom período no desenvolver deste trabalho. Aos meus companheiros de trabalho, pela contribuição no desenvolvimento deste trabalho em discussões, como também em muitas risadas e ótimas conversas que deixaram o dia-a-dia mais agradável. E a todos aqueles que contribuíram direta ou indiretamente para a realização deste trabalho, agradeço de coração. O meu, Muito obrigada! i SUMÁRIO Página RESUMO..................................................................................................................... ii ABSTRACT ................................................................................................................ iii LISTA DE FIGURAS .................................................................................................. iv LISTA DE TABELAS ................................................................................................... v 1 INTRODUÇÃO ......................................................................................................... 1 2 REVISÃO BIBLIOGRÁFICA .................................................................................... 3 2.1 Erosão Hídrica e o Escoamento superficial ....................................................... 3 2.2 Poluentes associados com a erosão hídrica ...................................................... 5 2.3 Modelos para estimar as perdas de solo e nutrientes no sedimento erodido .... 9 2.4 Estatística clássica e geoestatística na ciência do solo ................................... 12 2.5 Krigagem .......................................................................................................... 19 2.6 Simulação estocástica ..................................................................................... 20 3 MATERIAL E MÉTODO ......................................................................................... 22 3.1 Descrição do local de estudo .......................................................................... 22 3.2 Aquisição de dados para o desenvolvimento do estudo ................................. 23 3.3 Aquisição de dados para o WEPP ................................................................... 24 3.4 Classificação do sedimento erodido ................................................................. 24 3.5 Modelagem do fósforo transportado no sedimento erodido ............................. 25 3.6 Análises de desempenho dos modelos testados na estimativa de fósforo nos sedimentos erodidos .............................................................................................. 27 3.7 Análise Geostatística ....................................................................................... 28 4 RESULTADOS E DISCUSSÃO ............................................................................. 32 5 CONCLUSÃO ......................................................................................................... 46 REFERÊNCIAS ......................................................................................................... 47 ii MODELAGEM DA PERDA DE FÓSFORO POR EROSÃO HÍDRICA RESUMO- Com a crescente preocupação com as perdas de fósforo na água de enxurrada e sedimentos erodidos enriquecidos com P, que podem aumentar o risco de eutrofização dos corpos d´água superficiais, o interesse no uso de modelos que avaliam os impactos do uso do solo, tornou-se importante. O modelo WEPP, tem uma ampla gama de aplicabilidade, uma vez que pode ser usado para simular processos de erosão, escoamento e transporte de elementos orgânicos e químicos. Entretanto, no Brasil, não há casos em que um algoritmo tenha sido associado ao WEPP para modelar o transporte de P. Assim, o objetivo deste trabalho foi avaliar o desempenho de modelos para estimar as perdas de fósforo disponível, no sedimento erodido predito pelo WEPP, em uma pequena bacia hidrográfica. O presente trabalho foi realizado em área localizada no município de Tabapuã, noroeste do Estado de São Paulo. Os dados para alimentar o banco de dados do WEPP basearam-se nos componentes: climático, hidrológico e desenvolvimento vegetal. Os sedimentos erodidos foram classificados em: frações primárias, areia, silte, agregados pequenos e grandes. Com tais frações foi calculada a taxa de enriquecimento por fósforo com o WEPP (WER). Outros dois algoritmos foram utilizados para estimar o fósforo transportado com o sedimento erodido: ln(ER) = 2,682 – 0,278 ln(Sed) e) ln(ER) = 2 - 0,2 ln(Sed). Para a validação das estimativas de P nos sedimentos, utilizaram-se os critérios: o coeficiente de Nash-Sutcliffe (NS), o erro quadrático médio (RMSE), o coeficiente de massa residual (CRM) e o índice de concordância (d). Os modelos testados para estimar as taxas de enriquecimento do sedimento erodido e enxurrada por P são eficientes, quando em uso com as predições do Water Erosion Prediction Project (WEPP). As taxas de enriquecimento obtidas com o WEPP (WER), equação ln(ER) = 2,682 – 0,278 ln(Sed) e ln(ER) = 2 - 0,2 ln(Sed), respectivamente, apresentaram coeficiente de Nash-Sutcliffe (NS) próximo de 1. As perdas de P com o sedimento e enxurrada apresentaram moderado grau de dependência especial (GDE), enquanto que a erosão estimada pelo WEPP apresentou alto GDE. Perdas de fósforo solúvel com a enxurrada acima de 0,02 mg L-1, valor crítico para eutrofização, podem ser obtidas em 81% da área da bacia com uma probabilidade superior a 75%. Palavras-chave: Taxa de enriquecimento, Modelo WEPP, Variabilidade espacial iii MODELING THE LOSS OF PHOSPHORUS BY WATER EROSION ABSTRACT – With the growing concern over phosphorus losses with the runoff water and eroded sediment enriched with P that can increase the risk of eutrophication of surface water bodies, interest in the use of models that assess the impacts of land use, has if important. The WEPP model has a wide range of applicability, since it can be used to simulate erosion, drainage and transport and organic chemicals. However, in Brazil, there aren’t cases in which an algorithm has been associated with the WEPP to model the transport P. Thus, the objective with this study was to evaluate the performance of mode loss of phosphorus with eroded sediment in agricultural environment to estimate phosphorus losses available in the sediment eroded predicted the WEPP in a small catchment. This work was carried out in an area located in the municipality of Tabapuã, northwest of São Paulo. The sediment losses were estimated with the WEPP model. Two approaches (P- empirical) estimated ER according to an empirical relationship, and the other approach used WER calculated by WEPP (P-WEPP). Models to estimate enrichment rates of the sediment (Sed) and runoff with the Water Erosion Prediction Project (WEPP), WER - WEPP, an equation ln(ER) = 2.682 - 0.278 ln(Sed) and ln(ER) = 2 - 0.2 ln (Sed) present coefficient of Nash-Sutcliffe (NS) next 1. keywords: enrichment ratio, WEPP model, spatial variability iv LISTA DE FIGURAS Figura Página Figura 1 Exemplo de semivariograma experimental............................. 16 Figura 2 Representação esquemática da área em estudo ................ 23 Figura 3 Fluxograma das etapas de identificação dos modelos testados .................................................................................. 28 Figura 4 Valores medidos e estimados para fósforo disponível no sedimento erodido, com valores de taxas de enriquecimento (WER e ER) e obtidos com os modelos testados................... 35 Figura 5 Padrão espacial obtido por Krigagem ordinária da erosão estimada com o modelo WEPP............................................. 37 Figura 6 Distribuição espacial obtida por simulação sequencial gaussiana das perdas de fósforo no sedimento erodido ....... 42 Figura 7 Distribuição espacial obtida por simulação sequencial da probabilidade para a concentração de fósforo disponível >0,020 mgL-1........................................................................... 44 v LISTA DE TABELAS Tabela Página Tabela 1 Características granulométrica e mineralógica da área de eestudo........................................................................................... 22 Tabela 2 Valores da quantidade de sedimento erodido, concentrações de fósforo (P) no solo e no sedimento e taxas de enriquecimento por P para 17 encostas na área de estudo....... 33 Tabela 3 Desempenho de três métodos de estimativa da taxa de enriquecimento por fósforo do sedimento erodido...................... 34 Tabela 4 Fósforo no solo, taxas de enriquecimento, fósforo no sedimento e na água da enxurrada, e variáveis da erosão do solo em função da curvatura do relevo ...................................... 36 Tabela 5 Estatística descritiva das taxas de enriquecimento, fósforo no sedimento e na água da enxurrada, e variáveis da erosão do solo para 626 pontos nas encostas convexas e côncavas..................................................................................... 39 Tabela 6 Modelos e parâmetros dos semivariogramas simples ajustados aos dados................................................................................... 41 Tabela 7 Concentrações do sedimento erodido, fósforo no solo e na água da enxurrada em 17 encostas da bacia............................. 43 Tabela 8 Percentagem da área da bacia por intervalo de probabilidade para concentração de fósforo disponível >0,020 mgL-1 na enxurrada ................................................................................... 45 1 1 INTRODUÇÃO A preocupação relacionada aos problemas ambientais vem crescendo de maneira rápida e acentuada, com um objetivo especial: a busca pela sustentabilidade dos ecossistemas com o desenvolvimento econômico baseado na exploração dos recursos naturais. Mundialmente, a qualidade da água representa, talvez, o principal problema ambiental, pois as condições insuficientes e a disponibilidade para os múltiplos usos afetam diretamente a economia do mundo. As causas da degradação ambiental, ainda que conhecidas, dificilmente são evitáveis e, dentre elas, destaca-se a erosão hídrica do solo, com a consequente produção e transporte de sedimentos devido à ação do escoamento superficial. O fenômeno de erosão tem sido uma preocupação constante em todas as situações relativas à gestão do uso do solo e da água, tanto para compreender e prever as perdas de solo, o processo de transporte de sedimentos, bem como manter a qualidade da água nas bacias hidrográficas. Os sedimentos são considerados importante parâmetro de avaliação da qualidade da água, pois sua deposição reduz o armazenamento de água em lagos e reservatórios, além disso, afeta negativamente ecossistemas aquáticos por transportar nutrientes, pesticidas e agentes patogênicos via escoamento. Dentre os nutrientes, o fósforo (P) é um elemento essencial na formação e desenvolvimento de todos os seres vivos, inclusive na produção e rendimento das culturas agrícolas. No entanto, a perda de fósforo, por meio do escoamento superficial e nos sedimentos enriquecidos, acelera o processo de eutrofização dos corpos d’águas, o que compromete o ambiente aquático, afetando as comunidades que dependem e vivem nesses ambientes. Existem diversos métodos de pesquisas para estudar e avaliar as ocorrências e consequências do processo de erosão hídrica, e os modelos de predição são fundamentais, pois podem ser utilizados para estimativas de perdas de solos e nutrientes ligados aos sedimentos erodidos, bem como propiciam informações essenciais sobre as relações entre erosão e sistemas de manejo das culturas. 2 Um modelo de predição de erosão é representado por um conjunto de equações, baseado em informações como clima, solo, topografia, hidrologia e manejo do solo, capaz de simular processos de erosão do solo associados com encostas individuais ou pequenas bacias. Entre os modelos de predição o WEPP (Water Erosion Prediction Project) é amplamente usado por pesquisadores em regiões com clima temperado, condições bem diferentes do clima tropical. Por este motivo, torna-se de fundamental importância a verificação da aplicabilidade deste modelo para as condições edafoclimáticas brasileiras. Tal verificação deverá ser realizada com o propósito de auxiliar no planejamento agrícola, para fins de conservação do solo e da água, respeitando-se os limites estabelecidos pela legislação. No Brasil, até o presente, na maioria dos estudos sobre erosão hídrica, o modelo WEPP tem sido utilizado tradicionalmente para monitorar e representar a erosão espacialmente e temporalmente, com detalhamento dos processos de transporte e da deposição de sedimentos em encostas ou pequenas bacias hidrográficas. Neste contexto, o desafio científico e tecnológico deste estudo está relacionado à adaptação e aplicação do WEPP para modelar o transporte do fósforo. Portanto, a hipótese deste trabalho é que a taxa de enriquecimento dos sedimentos erodidos por fósforo pode ser obtida por meio de um algoritmo, o qual poderá ser utilizado com o modelo WEPP para estimar as perdas de P. Isto permitirá verificar possíveis impactos, que as mudanças no uso e manejo dos solos agrícolas poderão gerar no processo de eutrofização de corpos d´água pela presença de P. Os objetivos deste trabalho foram: 1) avaliar o desempenho de modelo WEPP para estimar as perdas de fósforo disponível no sedimento erodido; 2) avaliar o padrão de variabilidade espacial das perdas de fósforo com o sedimento erodido. 3 2 REVISÃO BIBLIOGRÁFICA 2.1 Erosão Hídrica e o Escoamento superficial A erosão hídrica é um dos principais problemas relativos ao uso dos solos no mundo e no Brasil, pois envolve o desprendimento, carreamento e deposição de suas partículas e que, simultaneamente, transportam elementos orgânicos e químicos adsorvidos para o ambiente aquático (PELLEGRINI et al., 2008). O excesso desses elementos afeta a estrutura e o funcionamento dos ecossistemas, comprometendo a fertilidade do solo e os efeitos na ordem física, financeira e social, pois existe um grupo preponderante de agricultores que não faz uso adequado de técnicas de manejo e/ou conservação do solo (OLIVEIRA et al., 2010). O processo de erosão hídrica é condicionado pelos fatores: intensidade da chuva, classe de solo, topografia, cobertura e manejo e práticas conservacionistas de suporte (HUDSON, 1977). Além disso, esse processo se divide em erosão em entressulcos e erosão em sulcos, de acordo com atributos distintos do seu fluxo superficial, os quais regem a mecânica do processo erosivo em cada uma dessas formas (MEYER et al., 1975). A erosão em entressulcos, conhecida também por erosão superficial ou laminar, caracteriza-se pela desagregação, por meio do impacto das gotas de chuva e transporte das partículas de menor diâmetro pelo escoamento laminar da água (FOSTER, 1982). Estas partículas, ao serem transportadas, de forma suave e uniforme, evidenciam a seletividade do processo erosivo, proporcionando maior adsorção de elementos minerais aos sedimentos finos. Enquanto na erosão em sulcos predominam partículas de tamanho maior, e a desagregação e o transporte ocorrem pelo escoamento da água nos sulcos (FOSTER, 1982; FREITAS et al., 2008 ; BERTOL et al., 2010). Embora a perda de solo no local pela erosão em entressulcos seja menor do que a erosão em sulcos, aquela é considerada um dos tipos de erosão que produz um efeito agressivo a camada superficial do solo, pois muitas vezes é difícil de ser observada. Durante os eventos pluviométricos intensos, o escoamento, que chega 4 aos corpos d´agua, transporta partículas de solo em suspensão, sedimentos, fertilizantes, nutrientes, matéria orgânica, agrotóxicos, entre outros, sendo ele considerado a principal fonte difusa de poluição de mananciais hídricos (SANTOS et al., 2010). O excesso destes poluentes (nutrientes, agrotóxicos, matéria orgânica, etc.) no meio aquático pode ocasionar desequilíbrio ambiental e, tanto o escoamento, quanto os sedimentos atuam como vetores naturais e importantes na transferência desses elementos para a biota aquática, além de acarretar problemas como assoreamento e eutrofização que comprometem os cursos d’água (O'GEEN et al., 2010 ; HU et al., 2013). De um modo geral, a concentração de poluentes apresenta-se de forma diversificada e dependente de fatores como o uso e manejo do solo, a cobertura vegetal, as declividades acentuadas, a textura e a compactação do solo e, também, da intensidade e frequência das precipitações, fatores estes que influenciam na composição química, física e biológica das águas (BERTOL et al., 2007 ; MELLO, 2009 ; MOURA et al., 2009). Neste contexto, ressalta-se a importância da cobertura vegetal do solo, pois a velocidade do escoamento superficial diminui com o aumento da sua porcentagem sobre o solo (COSTA et al., 2013). A presença total ou parcial da cobertura do solo, seja por meio de dosséis vegetativos ou resíduos vegetais, contribui para um controle mais efetivo da erosão hídrica e diminuição da velocidade do escoamento (AMARAL et al., 2013). Além de contribuírem para que partes dos nutrientes adsorvidos pelas plantas retornem ao solo e sejam aproveitados pelas culturas subsequentes (MALUF et al., 2015). Martins Filho et al. (2009) avaliaram a quantidade de palha de cana-de-açúcar que poderia ser removida do solo e, ainda, fornecer proteção significativa no controle das perdas de solo e nutrientes por erosão hídrica. Paula (2015) verificou a necessidade da cobertura mínima da superfície do solo de 42%, para que não haja taxa de enriquecimento (ER) do sedimento erodido por fósforo superior a 1,0. Em um trabalho posterior, Valim et. al. (2016), observou-se que a presença de resíduos vegetais de cana-de-açúcar sobre a superfície do solo aumenta o intervalo de tempo necessário para o início do escoamento superficial. 5 As aplicações periódicas e superficiais de altas doses de fertilizantes, em sistemas de cultivos de plantio direto (SPD), aumentam as quantidades de nutrientes nas camadas superficiais do solo. Ao serem aplicados na superfície do solo e quando não são absorvidos pelas plantas, consequentemente, os fertilizantes são carreados pelo escoamento superficial e/ou adsorvidos aos sedimentos, afetando significativamente os ecossistemas (BERTOL et al., 2007 ; MARTINS FILHO et al., 2009). Logo, a erosão hídrica dos solos e a consequente produção de sedimentos têm sido preocupação constante em todas as situações relativas à gestão do uso do solo e da água, tanto para compreender e prever os processos de transporte de sedimentos, como para manter a qualidade da água das bacias hidrográficas. 2.2 Poluentes associados com a erosão hídrica Nas últimas décadas, os sedimentos erodidos têm sido identificados como um potencial poluente, o qual pode afetar significativamente a ecologia aquática. Logo, os sedimentos podem representar uma séria ameaça à biodiversidade. Neste contexto, o monitoramento da qualidade dos sedimentos é considerado necessário para a proteção dos ambientes aquáticos. Pesquisadores como Silva et al. (2012) e Hortellani et al. (2013) utilizaram abordagens ligadas aos Valores Guias da Qualidade dos Sedimentos (VGQS), como orientadores da qualidade e da toxicidade dos sedimentos ligadas à estrutura das comunidades aquáticas, estratégia proposta pelo órgão ambiental canadense, “Canadian Council of Ministers of the Environment” (CCME). No Brasil, as resoluções de nº 357/ 2005 e de nº 430/ 2011 do Conselho Nacional do Meio Ambiente (CONAMA), apresentam parâmetros físicos, químicos e microbiológicos para o enquadramento dos corpos hídricos, o que não possibilita uma visão global dos efeitos gerados pelos sedimentos nos sistemas aquáticos, devido à ausência de padrões que permitam avaliar a qualidade dos sedimentos (BRASIL, 2005; BRASIL, 2011). Os sedimentos são formados por partículas sólidas de diferentes tamanhos, e representam um importante poluente para águas superficiais, devido ao seu papel de transferência e de destino de muitas substâncias, tais como: incluindo nutrientes, 6 metais pesados, pesticidas e outros contaminantes orgânicos (WALLING e COLLINS, 2008 ; BERTOL et al., 2010 ; SILVA et al., 2012). Diversos estudos ambientais, sobre a respeito do acúmulo de nutrientes nos sedimentos para avaliar a qualidade do corpo hídrico, têm sido realizados (FAGNANI et al. 2011; BORTOLUZZI et al. 2013; CHAKRABORTY et al., 2014 ; BHADHA, LANG, DAROUB, 2016). Ao entrarem nos ambientes aquáticos, os sedimentos possibilitam alterações na turbidez (penetração da luz), na temperatura da água e no oxigénio disponível, o que proporciona efeitos prejudiciais graves na biodiversidade. A magnitude do efeito do acúmulo de nutrientes nos corpos aquáticos depende da sua concentração, da duração da exposição, da composição química, do tamanho das partículas dos sólidos, da transferência, da dispersão e seu excesso. Na cadeia alimentar, a presença desses nutrientes, em grandes quantidades, é considerada poluente de alto potencial degradador, podendo causar danos ao ecossistema (BILOTTA; BRAZIER, 2008; YAKUTINA, 2011, SILVA; GASPARETTO, 2016). As concentrações de nutrientes no sedimento são altamente dinâmicas devido às transformações que ocorrem entre a fração líquida e a sólida (partículas) do solo, ou seja, os nutrientes da solução podem se ligar ao soluto de acordo com a capacidade de adsorção e de concentração das partículas existentes no solo (BOSCOV, 1997). Adsorção é um processo físico-químico pelo qual os íons ou moléculas são atraídos para as partículas sólidas (minerais e orgânicas) dotadas de cargas elétricas e com dimensões, em geral, inferiores a 1 µm. As partículas minerais, especialmente as frações de menor tamanho (argilas, silte e matéria orgânica), transportadas com os sedimentos erodidos têm áreas superficiais e densidade de cargas específicas elevadas, de modo a aumentar o potencial de adsorção de nutrientes, agroquímicos e metais pesados (VAN RAIJ, 2011). As partículas minerais possuem tamanhos diferentes e, geralmente, desenvolvem na superfície cargas negativas, que atraem e retêm fracamente vários cátions (ANDREOLI et al., 2005). Devido ao fato citado, as concentrações de 7 nutrientes podem ser mais elevadas nos sedimentos do que na fase líquida do solo (QUINTON et al., 2001). Dentre os elementos que podem estar adsorvidos nos sedimentos destaca-se o fósforo (P). O P é um elemento essencial para processos biológicos, bioquímicos dos seres vivos e está intimamente relacionado à produtividade e à diversidade em todos os ecossistemas terrestres e aquáticos (MANZONI ; PORPORATO, 2011). Entretanto, o baixo teor de P nos solos tropicais e a remoção pelas culturas e pelos processos erosivos dos solos intemperizados restringem a sua disponibilidade para as plantas e os animais. Isto se deve aos vários processos geoquímicos do solo, tais como solubilização, complexação, adsorção e precipitação que determinam a mobilidade e destino desse elemento (CHINTALA et al., 2014). As formas de fósforo no solo podem ser agrupadas em dois grandes grupos: fósforo orgânico (Porg) e fósforo inorgânico (Pin). O grupo do Porg está relacionado com a matéria orgânica presente no solo (NASH et al., 2014). O Pin pode ser separado em duas partes, o fósforo dos minerais primários e o fósforo adsorvido, formando diferentes compostos e com distintos graus de estabilidade química (NASH et al., 2014). Tais características classificam o elemento em: lábil e não lábil; diferenciando- o pelo grau de estabilidade ou solubilidade, conferindo um equilíbrio dinâmico entre as formas com diferentes efeitos: retenção, solubilização, disponibilidades, absorção pelas plantas ou organismos do solo (TURNER et.al., 2003; SANTOS et al., 2008; SHEN et al., 2011). Assim a fração lábil mantém um equilibro estável com P em solução. Já o P adsorvido (não lábil) entra lentamente em equilíbrio com P-solução, isto é, quando a planta retira o P da solução, parte do fósforo não lábil é liberada para repor e manter o equilíbrio. Entretanto, ao longo do tempo, o fósforo lábil pode se transformar em formas mais recalcitrantes, passando da superfície da partícula para o seu interior, resultando em maior estabilidade e menor possibilidade de dessorção do fosfato (SANTOS et al., 2008 ; VAN RAIJ, 2011). Nos processos erosivos hídricos, as perdas de P ocorrem pela sua adsorção nas superfícies dos sedimentos ou em suspensão no fluxo de água (BARBOSA et al., 2009), tornando-o indisponível as plantas (CESSA et. al., 2009). Essa adsorção 8 de íons sobre a superfície mineral dos solos tem ampla atuação sobre a disponibilidade e mobilidade de nutrientes e contaminantes (GEBLER et al., 2012). A retenção ou fixação dos nutrientes, em especial o P, nas superfícies dos minerais, varia de acordo com o material de origem, a intensidade de intemperismo, a composição da solução do solo, o manejo do solo, bem como com as características dos minerais, ou seja, a cristalinidade, a área superficial específica, e a concentração de radicais OH protonados – variações estas que interferem expressivamente na adsorção de nutrientes (INDA JUNIOR; KÄMPF, 2005 ; SILVA NETO et al., 2008; CAMPOS et al., 2016). Em solos tropicais e subtropicais, predominam minerais que influenciam características como a estrutura, a cor, a adsorção de nutrientes e importantes processos geoquímicos. O processo de adsorção do P ocorre fortemente nos solos intemperizados, com predomínio de argilas do grupo caulinita, argilomineriais 1:1, óxidos e hidróxidos de ferro (Fe) e de alumínio (Al) (VAN RAIJ, 1991; CHINTALA et al., 2014). Dentre os óxidos e hidróxidos de ferro, goethita e hematita são os minerais de maior ocorrência nos solos tropicais. Estes óxidos e hidróxidos se apresentam na forma de minerais cristalizados a materiais amorfos e, devido a sua abundância e diversidade em grau de cristalinidade, influenciam características físicas e químicas dos solos. Quanto aos teores de óxidos de Al, também apresentam correlação com a adsorção de P, podendo ser maior ou equivalente à dos óxidos de ferro (MESQUITA FILHO; TORRENT, 1993; FONTES; WEED, 1996). Outro aspecto a salientar é que os minerais goethita e hematita são importantes indicadores pedogenéticos, pois a constituição deles é assegurada pelas condições ambientais e por permanecerem por longo tempo no solo (KÄMPF ; CURI, 2000 ; IDA JUNIOR ; KÄMPF, 2005). Fink et al. (2014) verificaram a disponibilidade de fósforo, em solos tropicais e subtropicais, relacionada à adsorção deste elemento pelas superfícies minerais. Já Campos et al. (2016) avaliaram a capacidade máxima de adsorção de P e o índice de sorção de P para solos brasileiros, com diferentes atributos químicos, físicos e mineralógicos. 9 Cabe ressaltar que determinados sistemas de produção agrícola são considerados poluidores, devido à aplicação excessiva de nutrientes ao solo e sem considerar a variabilidade espacial dentro de uma determinada área de manejo. O aumento no uso destes sistemas de produção, para satisfazer as necessidades de uma crescente população humana, resulta na contaminação dos ecossistemas aquáticos com substâncias químicas, por meio da pulverização de agrotóxicos, da infiltração e do escoamento superficial (ARIAS-ESTÉVEZ et al., 2008 ; NUTTENS et al., 2016). 2.3 Modelos para estimar as perdas de solo e nutrientes no sedimento erodido A adequação de ferramentas de simulação das perdas de solo e de transporte de sedimentos, via escoamento em solos erodidos, tornou-se necessária para avaliações confiáveis de predições de impactos nas bacias hidrográficas. Segundo Quinton et. al. (2010), regiões agrícolas perdem, por ano, significativas quantidades de nutrientes (30 milhões de t ano–1 de N e 17,5 milhões de t ano–1 de P e carbono), comprometendo a fertilidade, a produtividade do solo, os ciclos biogeoquímicos, o clima e a qualidade d’água. Neste contexto e tendo em vista a diversidade de métodos de pesquisas para se estudar e avaliar as ocorrências e as consequências do processo de erosão hídrica, os modelos de predição são fundamentais, pois podem ser utilizados para estudos de perdas de solos e transporte de sedimentos em diferentes cenários de manejo do solo, devido à união de conhecimentos e de dados gerados pelos modelos físicos, sem a necessidade de realizar testes com custos elevados e gastos excessivos de tempo no campo. Com esses modelos os conservacionistas serão conduzidos à escolha de práticas de conservação que se encaixem nas necessidades da área específica (AKSOY ; KAVVAS, 2005; PIERI et al., 2007; DE MARIA, 2010). Um modelo de predição da erosão é representado por um conjunto de equações constituído por fatores como: clima, solo, topografia, manejo do solo, perda de solo e sedimentos (FLANAGAN et al. 2012). Para as zonas temperadas, os modelos de predição da erosão mais usados são a Equação Universal de Perdas de 10 Solo (Universal Soil Loss Equation - USLE) e a Equação Universal de Perdas de Solo Revisada (Revised Universal Soil Loss Equation – RUSLE 2, (VAN DIJK et al., 2002). Na zona tropical, modelos como o Projeto de Predição de Erosão Hídrica (Water Erosion Prediction Project - WEPP) têm proporcionado bons resultados na predição das taxas de erosão (AMORIM et al., 2010). O modelo WEPP, desenvolvido na década de 90, tem sido utilizado para estimar a perda e deposição de solo e alcançar uma ampla aplicabilidade em vários ambientes (FLANAGAN; NEARING, 1995; MAHMOODABADI ; CERDÀ, 2013). Na literatura brasileira há registros de trabalhos com a aplicação do modelo WEEP; dentre os mais recentes, cita-se os de Amorim et al. (2010); Costa (2015); Ibiapina (2015) e Moraes (2016), realizados para algumas condições de manejo do solo. Nos referidos trabalhos, realizados no Brasil, apenas o desempenho do modelo WEPP e a sua aplicação para estimativas das perdas de solo e o transporte de sedimentos foram avaliados, para pequenas bacias hidrográficas. Portanto, não tiveram como objetivo estimar a perda de nutrientes adsorvidos aos sedimentos. No entanto, o modelo WEPP, em alguns países, tem sido testado para simular a perda de pesticidas e nutrientes, em especial as perdas de fósforo (PEREZ-BIDEGAIN et al., 2010; SAVABI et al., 2011). É sabido que a maioria dos nutrientes são transportados adsorvidos aos sedimentos erodidos (SILVA et al., 2012) e, ainda, que o WEPP tem capacidade para simular o movimento das partículas do solo e calcular as taxas de enriquecimento em áreas com até 260 ha (PEREZ-BIDEGAIN et al. 2010 ; WANG et al., 2010). As perdas de fósforo com a água da enxurrada e sedimentos erodidos enriquecidos com P podem aumentar o risco de eutrofização dos corpos d´água superficiais (WANG et al., 2011). Em muitos casos, a eutrofização restringe o uso da água para a pesca, recreação, uso industrial e produção agrícola, devido à proliferação de algas e plantas aquáticas (WANG et al., 2011; SINGH et al., 2011). Deste modo, é crescente a preocupação em reduzir a quantidade de P transferida para os corpos d´água. Neste contexto, o interesse no uso de modelos que avaliam os impactos de diferentes práticas de manejo do solo no transporte de fósforo com a erosão torna-se importante. 11 Autores como Sharpley e Kleinman (2003) realizaram uma pesquisa, na qual estabeleceram uma relação logarítmica entre o grau de enriquecimento por fósforo com a quantidade de sedimentos erodidos. Tal relação foi estabelecida para bacias hidrográficas apresentando solos com diferentes texturas, declividades e usos com culturas. Diante do exposto, a combinação de um algoritmo de transporte de fósforo com o WEPP poderá expandir a sua aplicabilidade na avaliação da qualidade da água. Além disso, poderá permitir investigar impactos ambientais em função dos usos e manejos dos solos nas bacias hidrográficas. Devido a este fato, um algoritmo foi proposto por FRERE et al. (1980), para estimar a quantidade de P transportado com os sedimentos erodidos de pequenas áreas ou bacias hidrográficas como: PSed = PSolo ₓ ER ₓ Sed (1) em que, PSed é quantidade de P no sedimento (kg ha-1), PSolo é o conteúdo de fósforo no solo (kg kg-1), ER é a taxa de enriquecimento por fósforo e Sed é a quantidade de sedimento erodido (kg ha-1). Um inconveniente da Equação 1 é que a ER necessita ser estimada, ao contrário do conteúdo de fósforo no solo e da quantidade de sedimento erodido, que podem ser medidos. Deste modo, foi proposto por Menzel, (1980) modelar a ER como: ln(ER) = a + b ₓ ln(Sed) (2) em que, a e b são parâmetros empíricos. No WEPP, o sedimento derivado de encostas ou pequenas bacias hidrográficas é predito em cinco classes de tamanho de partículas: frações primárias argila, silte e areia; agregados pequenos constituídos por argila, silte e matéria orgânica; agregados grandes constituídos pelas três frações primárias e matéria orgânica (ELLIOT et al., 2015). A fração de sedimentos erodidos, em cada uma das cinco classes de tamanho, derivados a partir de uma encosta ou de um divisor de águas é, portanto, gerada pelo modelo WEPP. O WEPP calcula a taxa de enriquecimento (WER) com base na superfície específica (SS), a qual é determinada como o valor da SS do sedimento erodido dividida pelo valor da SS do solo originalmente presente na encosta. A WER pode ser utilizada para determinar o aumento da concentração de um elemento, poluente 12 ou não, presente no sedimento em comparação com a sua concentração no solo original de uma encosta, tal como demonstrado nos trabalhos de Elliot et al. (2015); Wilson et al. (2016) e Papanicolaou et al. (2015). A WER pode ser determinada como: (3) em que, WER é a taxa de enriquecimento, SSASed é a superfície específica do sedimento (m2 g-1); SAASolo é a superfície específica do solo original (m2 g-1). A taxa de enriquecimento, computada no WEPP, pode ser testada para estimar a quantidade de fósforo transportado com o sedimento como: PSed = PSolo ₓ WER ₓ Sed (4) em que, PSed é quantidade de P no sedimento (kg ha-1), PSolo é o conteúdo de fósforo no solo (kg kg-1), WER é a taxa de enriquecimento por fósforo, computada pelo WEPP e Sed é a quantidade de sedimento erodido (kg ha-1). No Brasil, o modelo WEPP tem sido avaliado para monitorar e representar a erosão espacialmente e temporalmente, com detalhamento dos processos de transporte e da deposição de sedimentos, em encostas ou pequenas bacias hidrográficas. Contudo, não há estudos em que um algoritmo tenha sido associado ao WEPP para modelar o transporte de P. Desde 1995 há um interesse em adicionar ao WEPP uma rotina ou algoritmo para o transporte de substâncias químicas (PEREZ-BIDEGAIN et al., 2010), contudo isto ainda não ocorreu até o presente. 2.4 Estatística clássica e geoestatística na ciência do solo A erosão do solo é um fenômeno dinâmico, sujeito às variações temporais e espaciais. Predições acuradas do processo de erosão requerem suficiente entendimento da variabilidade de fatores como: erosividade das chuvas, erodibilidade do solo, relevo, cobertura vegetal, manejo do solo e práticas de suporte. Portanto, em questões ambientais relativas às perdas de substâncias químicas por processos erosivos do solo, é relevante considerar a variabilidade da erosão em diferentes escalas espaciais e temporais (MARTINS FILHO, 2007). 13 Na área da Ciência do Solo, tem sido muito comum a utilização da estatística clássica no exame de resultados experimentais, em três classes: 1) medidas de posição ou tendência central; 2) medidas de dispersão; e 3) medidas da forma da distribuição. Dentre as medidas de tendência ou posição central, destacam-se a média e a mediana, sendo a média aritmética dos dados a mais utilizada. Medidas de tendência central são, normalmente, insuficientes para descrever completamente um conjunto de dados (MARTINS FILHO, 2007). Faz-se necessário, portanto, o uso de medidas de dispersão como desvio-padrão, coeficiente de variação (CV), valores máximo e mínimo. O desvio-padrão é igual à raiz quadrada da variância, o que corresponde à média das diferenças quadráticas dos valores observados em relação a sua média. Portanto, caso os valores tendam a se concentrar próximos da média, a variância é baixa e, em caso contrário, é alta confirmando observações de Warrick e Nielsen (1980). O CV é determinado dividindo-se o desvio-padrão pela média, o qual indica o grau de variabilidade da variável em análise. Contudo, o CV não reflete o significado físico da variável. Em geral, os atributos e a erosão do solo têm alto CV (NEARING et al., 1999; ELLIOT et al., 2001). Na avaliação, se os dados experimentais seguem uma distribuição normal, pode-se utilizar os coeficientes de assimetria e curtose. Idealmente devem ter valor zero ou o mais próximo disso. Valores entre -2 e +2 para assimetria e curtose podem ser aceitos, segundo Ortiz (2003). Assimetria e curtose são medidas muito utilizadas em sedimentologia, mas também são guias úteis no momento de se definir pelo tipo de distribuição estatística dos dados de estudo. As propriedades do solo, em geral, não apresentam distribuição normal. Elas não estão aleatoriamente distribuídas no espaço, mas obedecem a um arranjo estrutural (MARTINS FILHO, 2007). Tal arranjo tem uma dimensão característica que é o seu domínio, que corresponde à distância dentro da qual há interdependência dos valores observados (MARTINS FILHO, 2007). Quando a distribuição dos dados não é normal, a média aritmética sofre influência dos valores extremos e esta não representará adequadamente o conjunto de dados. Neste caso, a estatística clássica pode não ser a ferramenta mais adequada. 14 Para utilizar a estatística clássica é necessário considerar que o solo seja homogêneo. Há, ainda, o problema da impossibilidade de se saber, previamente à amostragem, se as amostras são dependentes ou independentes. Como os solos são heterogêneos, é necessária a utilização de métodos estatísticos adicionais que considerem tais problemas. Desse modo, para a análise da maioria das variáveis, oriundas de fenômenos naturais, é necessário se admitir alguma incerteza no comportamento destes fenômenos entre duas posições espaciais da amostragem (ISAAKS ; SRIVASTAVA, 1989). Dentro deste contexto, a geoestatística surge como uma alternativa consistente para modelar este comportamento, por meio do uso de funções aleatórias. O uso da geoestatística pode ser uma ferramenta fundamental para as predições da erosão e tomadas de decisões em atividades de planejamento, visando à conservação do solo e da água (MARTINS FILHO, 2007). A geoestatística pode ser descrita como um conjunto de técnicas que permitem modelar a variação no espaço e/ou no tempo de uma série de dados, a qual é também definida como a teoria das variáveis regionalizadas (JOURNEL ; HUIJBREGTS, 1991). Uma variável regionalizada, Z(xi), é uma função espacial numérica, que varia de um local para outro, com continuidade aparente e cuja variação não pode ser representada por uma expressão matemática simples. Essa variável Z(xi) assume diferentes valores de Z dada a localização ou a posição x numa área (CRESSIE, 1991). A base da geoestatística está centrada nos conceitos de variáveis regionalizadas, funções aleatórias e estacionariedade. Uma variável aleatória é uma particular medida do que se espera variar de acordo com uma lei probabilística, a qual é caracterizada por parâmetros da distribuição como valor esperado e variância (TRANGMAR et al., 1985). Quando uma variável aleatória assume valores diferentes em função da posição de amostragem no campo, define-se uma variável regionalizada. O conjunto de todas as possíveis realizações da variável aleatória, em todas as posições de uma área, constitui a função aleatória. A continuidade da função aleatória é definida como estacionariedade. Assume-se, portanto, que as propriedades estatísticas (média, variância, covariância, momentos de maior ordem, etc.) são estacionárias no espaço, isto é, 15 não variam com a translação (COSTA, 1999). Trangmar et al. (1985) definiram a hipótese de estacionariedade de primeira ordem como sendo aquela em que o valor esperado da função aleatória Z(xi) é o mesmo para toda a área, independentemente da posição que ocupa ou da distância de separação (h), tal que: E[Z(xi)] = m (5) E[Z(xi) – Z(xi+h)] = 0 (6) A estacionariedade de segunda ordem, também conhecida como estacionariedade forte, é admitida quando a função aleatória atende à estacionariedade de primeira ordem e à co-variância espacial C(h), para cada par de valores Z(xi), Z(xi+h) separados por uma distância h, é igual em toda a área estudada e depende apenas de h, o que implica a ocorrência de variância finita (COSTA, 1999). Esta autora considera que a co-variância C(h) de cada Z(xi) e Z(xi+h) é a mesma, em qualquer posição da área de estudo, tal que: C(h) = E[Z(xi)-m] [Z(xi+h)] = 0 (7) Uma forma de estacionariedade menos restritiva, a denominada hipótese intrínseca, pode ser assumida (SOUZA, 2004). Esta hipótese requer que a variância Z(xi) – Z(xi+h) seja finita, para todos os vetores h, e independentemente da posição na área em análise, tal que: VAR[Z(xi) – Z(xi+h)] = E[Z(xi) – Z(xi+h)]2 = 2(h) (8) A semivariância é a medida do grau de dependência espacial entre duas amostras, a qual é representada por (h). A semivariância depende do vetor de separação h. Idealmente (h) é zero quando h=0 (TRANGMAR et al., 1985). Na análise geoestatística, a variabilidade espacial é profundamente avaliada e modelada, para em seguida se empregar técnicas apropriadas de estimativas, cujos resultados serão imagens representativas da distribuição no espaço, das propriedades que estão sendo analisadas (STURARO, 1993). A geoestatística estuda o grau de dependência no espaço da grandeza medida e o alcance ou o domínio de cada amostragem, considerando a posição relativa da amostra no terreno ou sua distribuição espacial, o que a difere da estatística clássica que, devido ao uso das médias, das estimativas da variância e do coeficiente de variação, não o faz (JOURNEL ; HUIJBREGTS, 1991; ISAAKS ; SRIVASTAVA, 1989). 16 A utilização de técnicas geoestatísticas permite a elaboração de semivariogramas, que estabelecem a dependência espacial de parâmetros do solo, possibilitando a construção de estimativas de médias e o dimensionamento de amostras com erros mínimos (GUIMARÃES et al., 1992). Segundo Vieira et al. (1983), o semivariograma é um gráfico que caracteriza a estrutura de dependência espacial da(s) variável(eis) sob estudo, ou seja, é uma função que relaciona a semivariância com o vetor distância, podendo ser representado analítica ou graficamente. A semivariância pode ser estimada pela equação: (h) = N(h)2 1      N(h) 1i 2 ii hxZ)(x Z (9) em que,  (h) – semivariância estimada; N (h) – número de pares de valores observados Z (xi), Z (xi + h) separados pela distância h. O gráfico do variograma experimental é caracterizado por três parâmetros, representados na Figura 1: a) O efeito pepita (C0); b) O patamar (C0 + C1); c) O alcance (a). Figura 1. Exemplo de semivariograma experimental, com os parâmetros C0 (efeito pepita), C0+C1 (patamar) e (a) alcance. O chamado efeito pepita (C0) é a descontinuidade do semivariograma próximo à sua origem, e expressa tanto a variabilidade para uma escala menor que o intervalo de amostragem quanto o erro de mensurações. O efeito pepita, que V a r iâ n c ia e s t ru tu ra l ( C 1 ) P a ta m a r ( C o + C 1 ) A lc a n c e (a ) Va ri â n cia p e p i ta ( C o ) S e m iv a r iâ n c i a M o d e lo a ju s t a d o D is tâ n c ia (m ) S em iv ar iâ nc ia 17 representa a variabilidade não-explicada, podendo ser originada dos erros de medição ou de microvariações não-detectadas quando considerada a distância de amostragem utilizada (CAMBARDELLA et al., 1994; SALVIANO et al., 1998; OLIVEIRA et al., 1999). Quando os dados apresentam uma variabilidade aleatória, que não se ajusta a nenhum modelo matemático, diz-se que ocorre o efeito pepita puro. Para distâncias superiores ao alcance, o patamar (Co+C1), quando existe, ocorre a partir da distância na qual a variância aproxima-se assintoticamente do máximo. O alcance (a) expressa a distância além da qual as amostras não estão correlacionadas (TRANGMAR et al., 1985). A razão entre o efeito pepita (C0) e o patamar (Co+C1) tem sido utilizada para definir o grau de dependência espacial de uma variável, segundo critérios de Cambardella et al. (1994), como: 1) Co/(Co+C1)  25%, variável com forte dependência espacial; 2) Co/(Co+C1) entre 25% e 75%, variável com moderada dependência espacial; 3) Co/(Co+C1) > 75%, variável com fraca dependência espacial. Os semivariogramas apresentam formas variáveis em função dos dados e do intervalo amostral utilizado (OLIVER e WEBSTER, 2014). Cada variável analisada é submetida ao ajuste de modelos matemáticos que irão definir os semivariogramas. O modelo matemático deve ser ajustado o mais próximo possível dos dados observados (McBRATNEY ; WEBSTER, 1986). Para Vieira (2000), o ideal é que a semivariância aumente com a distância entre os pontos amostrados até atingir um valor mais ou menos constante. Uma maior dependência espacial e um melhor ajuste são conseguidos quando, no processo de modelagem, o índice de melhor ajuste (IMA) é próximo de zero (PANNATIER, 1996). Conforme McBratney e Webster (1986), os semivariogramas são ajustados a partir de padrões gráficos preestabelecidos, como segue: a) Modelo esférico: (h) = C0 + C1 3/2 (h/a) – 1/2 (h/a)3  0 < h < a (10) (h) = C0 + C1 h  a (11) 18 Este modelo é obtido selecionando-se os valores do efeito pepita e do patamar e, depois, passando-se uma reta que intercepte o eixo y em C0 e seja tangente aos primeiros pontos próximos de h = 0. Tal reta cruzará o patamar a uma distância a’= (2/3)a. Logo, o alcance será a = 3a’/2. Segundo Vieira (2000), o modelo esférico é linear até (1/3)a. Este é o modelo que mais se ajusta aos atributos do solo (TRANGMAR et al., 1985). b) Modelo exponencial: (h) = C0 + C1 1- exp (-3 h/a) 0 < h < d (12) Já no modelo exponencial, o patamar é atingido assintoticamente, com o alcance prático definido como a distância na qual o valor do modelo é 95% do patamar (ISAAKS ; SRIVASTAVA, 1989). Esta maneira de atingir o patamar é o que difere o modelo exponencial do modelo esférico. Segundo Vieira (2000), os parâmetros C0 e C1 são determinados da mesma maneira que para o modelo esférico. c) Modelo gaussiano: (h) = C0 + C1 1- exp (-3 (h/a)2) 0 < h < d (13) O modelo gaussiano é usado, muitas vezes, para modelar fenômenos extremamente contínuos (ISAAKS ; SRIVASTAVA, 1989). O que caracteriza este modelo é o seu ponto de inflexão próximo à origem. Ele também atinge o patamar assintoticamente. O parâmetro a é definido como o alcance prático ou distância na qual o valor do modelo atinge 95% do patamar (ISAAKS ; SRIVASTAVA, 1989). d) Modelo linear (h) = C0 + (C1/a) h 0 < h < 0 (14) (h) = C0 + C1 h  a (15) em que, C1/a é o coeficiente angular para 0 < h < a. O coeficiente angular é determinado pela inclinação da reta que passa pelos primeiros pontos de (h), dando-se maior peso àqueles que correspondem ao maior número de pares. Neste modelo, o patamar é determinado por inspeção. O efeito pepita, C0, é determinado 19 pela interseção da reta no eixo  (h). Já o alcance (a) é o valor de h correspondente ao cruzamento da reta inicial com o patamar e C1 = patamar - C0 (VIEIRA, 2000). Escolher o modelo mais adequado não é um procedimento automático conforme apontaram Oliver e Webster (2014). Em geoestatística, é comum o ajuste visual do modelo selecionado aos pontos experimentais, o que carece de sustentação estatística. Segundo McBratney e Webster (1986), ajustes são normalmente feitos pelos métodos da minimização do quadrado dos desvios, assumindo normalidade e independência dos resíduos e homogeneidade de variâncias. O uso desse método para ajustar semivariogramas é muito criticado em função dessas hipóteses assumidas, sobretudo no que diz respeito à variância. O método dos mínimos quadrados generalizados contorna essa questão (CRESSIE, 1991). Uma alternativa tem sido a utilização do método dos mínimos quadrados ponderados, com o qual o peso associado a cada ponto experimental é uma função do número de pares de observações que o originou. 2.5 Krigagem Numa análise geoestatística, nem sempre o interesse é apenas a obtenção de um modelo de dependência espacial, mas também é a predição de valores da variável em apreço para locais ou pontos não amostrados numa área. O interesse pode estar centrado em um ou mais pontos específicos da área ou numa malha de pontos interpolados que permitam visualizar o comportamento da variável numa região através de um mapa de isolinhas ou de superfície (ORTIZ, 2006). Para tanto, é necessária a aplicação de um método de interpolação. Um método de interpolação bastante utilizado é a krigagem, a qual é uma técnica desenvolvida para estimar valores para locais não-amostrados, a partir da correlação espacial definida pelo semivariograma. Segundo Isaaks e Srivastava (1989), é um modelo de interpolação utilizado para estimar valores de uma propriedade em locais não-amostrados, a partir de valores vizinhos resultantes da amostragem realizada. Esse interpolador pondera os vizinhos do ponto a ser estimado, segundo critérios de não-tedenciosidade e variância mínima (GONÇALVES, 1997). Segundo 20 LI e HEAP (2008), existem vários tipos de krigagem: simples, ordinária, universal, indicadora, probabilística etc. Com a interpolação de valores para cada local não amostrado, é possível estabelecer, pela krigagem, um mapeamento da área em estudo para as propriedades dos solos sob investigação, inclusive com a construção de linhas de isovalores (BURGESS; WEBSTER, 1980). A possibilidade de utilização dessas ferramentas da geoestatística representa um ganho científico muito importante, em especial para a ciência do solo, como um elemento auxiliar no planejamento das atividades agrícolas e no fornecimento de subsídios para tomadas de decisões estratégicas visando ao uso, ao manejo e à conservação dos recursos naturais solo e água. 2.6 Simulação estocástica A simulação estocástica é um processo de gerar múltiplas e equiprováveis realizações de variáveis aleatórias componentes de um modelo de uma função aleatória (DEUTSCH; JOURNEL, 1998). As realizações {z(l)(u), u  A, com l=1,…,L}, geralmente dispostas segundo uma malha regular, representam as L possibilidades da distribuição espacial dos valores do atributo z(u) ao longo de uma área A, tal que cada realização é denotada com o índice (l) . Cada uma das realizações ou “imagens estocásticas” deve refletir as propriedades impostas no modelo da função aleatória z(u). Portanto, quanto mais propriedades forem inferidas dos dados amostrais e incorporadas no modelo de função aleatória z(u), tanto melhor será este modelo. As variáveis aleatórias (VA) de um campo aleatório, que representa a distribuição de um atributo espacial, não são independentes (FELGUEIRAS, 1999). Há uma correlação entre essas variáveis que, sob a hipótese de estacionariedade de segunda ordem, depende do vetor de separação, h, determinado pela posição espacial das amostras. Deste modo, a simulação conjunta de m variáveis aleatórias, que compõem um campo aleatório, requer a inferência e a modelagem da matriz de covariância cruzada entre essas variáveis. Entretanto, este procedimento é muito custoso computacionalmente e, por isso, não é adotado na prática (FELGUEIRAS, 1999). Deutsch e Journel (1998) apresentaram uma metodologia eficiente que 21 aproxima a simulação conjunta de muitas variáveis aleatórias, denominada simulação sequencial condicionada. Essa simulação usa a função de distribuição acumulada condicionada (fdac) às VA mais correlacionadas, para obter valores z(u) de uma variável aleatória z em cada posição u  Α. Segundo Chilès e Delfiner (1999), uma função aleatória tem um número infinito de realizações, entre as quais as simulações condicionais, que assumem nos pontos amostrais os mesmos valores que aqueles observados. O condicionamento considera os dados amostrais originais e também os valores pré-simulados dentro da vizinhança u. Esse condicionamento confere certa robustez às simulações (JOURNEL; HUIJBREGTS, 1991), que podem assim ser consideradas como possíveis representações da variável regionalizada z(u). As realizações equiprováveis zs(u) diferirão ponto a ponto (ou bloco a bloco) da realização real z(u), mas, estatisticamente no conjunto global, mostrarão a mesma estrutura de variabilidade ou mesma estatística espacial, i.e., mesmos histogramas e variogramas (JOURNEL; HUIJBREGTS, 1991). Essas diferentes realizações podem ser vistas como possíveis cenários do mesmo fenômeno caracterizado pela função aleatória z(u). Adicionalmente, com as realizações simuladas há a vantagem de serem conhecidos todos os pontos e não apenas os pontos amostrais. A série de realizações fornece uma medida visual e quantitativa, ou seja, um modelo, de incerteza espacial (GOOVAERTS, 1997). 22 3 MATERIAL E MÉTODO 3.1 Descrição do local de estudo A área utilizada neste trabalho localiza-se no município de Tabapuã, noroeste do Estado de São Paulo, cujas coordenadas geográficas são: Latitude 21º 05’ 57 S e Longitude 49º 01’ 02 W (Figura 2a). O clima da região, segundo a classificação de Thornthwaite, é do tipo megatérmico, C2dA’a’, subúmido chuvoso com pequeno ou nenhum excedente de água, e evapotranspiração de verão menor que 48% do total anual (THORNTHWAITE, 1948). A precipitação média anual é de 1.318 mm. A vegetação primária da região foi classificada como floresta pluvial estacional e cerrado. A área tem 220 ha e apresenta um histórico de mais de 20 anos consecutivos com cultivo de cana-de-açúcar, sendo mais de 10 anos no sistema de colheita mecanizada sem queima. Logo, a cobertura vegetal da área, na época da realização deste trabalho, era constituída por resíduos de cana-de-açúcar. O solo da área foi classificado por Sanchez et al. (2009) como Argissolo Vermelho eutrófico, conforme os critérios da Embrapa (2013). O material de origem foi identificado como rochas areníticas sedimentares do Grupo Bauru (IPT, 1981). Foram identificadas duas formas de relevo, uma côncava e outra convexa (Figura 2b). A composição granulométrica e mineralógica estão representadas na Tabela 1, de acordo com Barbieri et al., (2009); Camargo et al. (2010); Silva Junior et al. (2012); Camargo et al.(2013a). Tabela 1 – Características granulométrica e mineralógica da área de estudo. Relevo Argila Silte AT MO Gt Ht Gb/(Ct+Gb) Fe2O3 ----------------------- g kg-1 ---------------------------- ---%--- Côncava 171 35 796 1,6 10,52 12,39 0,12 2 Convexa 190 43 784 1,5 13,24 23,20 0,06 4 AT- Areia total; MO- Matéria Orgânica; Gt-Goethita; Gb-Gbisita; Ct-Caulinita; Ht–Hematita; Fe2O3 - Ferro total 23 3.2 Aquisição de dados para o desenvolvimento do estudo Os dados medidos para a estimativa e validação dos modelos foram obtidos de Martins Filho (2007), que estimou a erosão para 626 pontos, espaçados de 50 em 50 m com o WEPP. Das 89 encostas selecionadas por Ibiapina (2015), selecionou-se 17 pontos de observação (17 encostas), dos quais obteve-se valores de concentração de fósforo no solo (Psolo) e no sedimentos erodidos (Psed). Assim, com a concentração de Psed sobre Psolo obteve-se a taxa de enriquecimento (ERsed) observada. As ERsed serviram de referência para validar os modelos testados. Figura 2. Localização da área próximo ao Rio Turvo (a), mapa de declividade e 17 pontos de observação para coleta de informações da chuva simulada (b), simulação do fluxo superficial de água (c) e mapa planialtimétrico com as cotas de altitude em metros (d). Adaptado de Ibiapina (2015). 24 3.3 Aquisição de dados para o WEPP Foram utilizadas 5 tipos de informações de entrada no modelo WEPP: (i) climática, (ii) hidrológica, (iii) fisiologia das culturas de cana-de-açúcar, (iv) tipo de solo e (v) dados de erosão/deposição, (Figura 2) (CHAVES, 1992). A informação climática foi feita estocasticamente, e utilizado modelo paralelo CLIGEN para gerar dados a partir de observações de precipitação, temperatura, ponto de orvalho e radiação solar da estação meteorológica do município de Catanduva – SP, onde os dados são mais completos e correspondentes à realidade de Tabapuã. A informação hidrológica foi gerada a partir das características da topografia, infiltração e evapotranspiração (Figura 2b, 2c e 2d). As informações do solo (textura, densidade e conteúdo de matéria orgânica) são do banco de dados de Martins Filho (2007); Sanchez et al. (2009). A classificação granulométrica nas informações de entrada sobre o solo seguiram os critérios de classificação do USDA (2014). As informações sobre características fisiológicas das culturas de cana-de açúcar foram extraídas do pacote utilizado pelo modelo WEEP. Nos cenários do WEPP, as coberturas do solo foram equivalentes a 90% e para as informações de entrada sobre o manejo foram utilizados os valores descritos por Pieri et al. (2007). 3.4 Classificação do sedimento erodido A composição do sedimento erodido foi estabelecida em cinco classes: 1) frações primárias (<0,002 mm), areia (0,06 - 2 mm), e silte (0,002-0,06 mm); 2) agregados pequenos constituídos por argila, silte e matéria orgânica; agregados grandes constituídos pelas três frações primárias e matéria orgânica (ELLIOT et al., 2015). Para o fracionamento do material erodido, a abordagem proposta por Foster et al. (1985) e incorporada no WEPP foi a utilizada. Tal fracionamento permitiu obter a taxa de enriquecimento WER para testes de validação em 17 encostas, das 89 estabelecidas na área de estudo (Figura 2b). 25 A fração argila primária nos sedimentos foi obtida como: Fcl = 0,26 Ocl (16) em que: Fcl é a fração de argila primária no sedimento e Ocl é a fração argila no solo original. As frações dos agregados pequenos e silte primário foram estimadas como: Fsg = 1,8 Ocl Ocl < 0,25 (17) Fsg = 0,45 – 0,6 (Ocl – 0,25) 0,25 ≤ Ocl ≤ 0,50 (18) Fsg = 0,6 Ocl > 0,50 (19) Fsi = Osi – Fsg (20) em que, Fsg é a fração agregados pequenos, Fsi é a fração silte primário, Osi é a fração silte no solo original. Já a fração areia primária foi calculada como: Fsa = Osa (1 – Ocl) 5 (21) em que, Osa = fração de areia no solo original. Uma vez que a soma das frações para todas as classes de partículas deve ser igual a um, a fração da classe dos agregados grandes foi calculada como: Flg = 1 – Fcl – Fsi – Fsg – Fsa (22) em que, Flg é a fração agregados grandes. Se Flg, a partir da Equação (22), é negativo, valores para as outras frações são proporcionalmente reduzidos para F lg dar zero. 3.5 Modelagem do fósforo transportado no sedimento erodido A estimativa do fósforo transportado com os sedimentos erodidos foi obtida adotando-se dois algoritmos. No primeiro algoritmo utilizou-se: PSed = PSolo ₓ ER ₓ Sed (23) Ln(ER) = a + b ₓ Ln(Sed) (24) em que, PSed é quantidade de P no sedimento (kg ha-1), PSolo é o conteúdo de fósforo no solo (kg kg-1), ER é a taxa de enriquecimento por fósforo e Sed é a quantidade de sedimento erodido (kg ha-1). Os valores inicialmente utilizados de a e b foram 2,0 e - 0,2 (MENZEL, 1980). Também os valores de a e b foram calibrados com valores de 26 ER obtidos em 7 encostas e, posteriormente, a equação resultante teve o seu desempenho testado utilizando-se a ER de 10 encostas. Os valores de ER medidos foram obtidos dividindo-se a concentração de P no sedimento erodido pela sua concentração no solo original. Para tanto, utilizou-se resultados experimentais, não publicados, obtidos por Martins Filho (2007). No segundo algoritmo, a quantidade de fósforo transportado com o sedimento erodido foi estimada como: PSed = PSolo ₓ WER ₓ Sed (25) WER (26) em que, PSed é quantidade de P no sedimento (kg ha-1), PSolo é o conteúdo de fósforo no solo (kg kg-1), WER é a taxa de enriquecimento por fósforo computada pelo WEPP e Sed é a quantidade de sedimento erodido (kg ha-1), em que, WER é a taxa de enriquecimento, SSASed é a superfície específica do sedimento (m2 g-1); SAASolo é a superfície específica do solo original (m2 g-1). O fluxograma das etapas de identificação dos modelos testados está representado na Figura 3. 27 Figura 3 – Fluxograma das etapas de identificação dos modelos testados 3.6 Análises de desempenho dos modelos testados na estimativa de fósforo nos sedimentos erodidos O desempenho dos modelos estudados foi avaliado utilizando-se diferentes critérios. Os critérios incluíram o coeficiente de Nash-Sutcliffe (NS), o erro quadrático médio (RMSE), valores das médias dos erros ou diferenças médias (MD) e o índice de concordância (d), determinados como: NS = 1 (27) RMSE (28) 28 MD = )O(P i n 1i i   / n (29) d = 1 (30) em que, Pi são valores estimados pelo modelo, Oi são os valores observados ou medidos, n é o número de valores observados e O a média dos valores medidos, respectivamente O valor máximo para o coeficiente NS é 1. O MD e o NS podem ser negativos. Se o coeficiente de NS for menor que 0 (zero) o modelo prediz valores piores do que o simples uso da média observada. Quando o MD apresentou sinal (+) ou (-), foi assumido que os valores preditos, em média, superestimavam ou subestimavam os valores observados, respectivamente. O índice d é delimitado entre 0 e 1. 3.7 Análise Geostatística A análise geoestatística foi realizada seguindo-se a metodologia descrita em Vieira et al. (1983) e Robertson (1998), em três etapas: 1) análise espacial para obter o semivariograma; 2) escolha de modelo para ajustar o semivariograma e 3) obtenção de valores das variáveis em estudo, em local não-amostrado, para interpolação através da krigagem. Para ajuste do semivariograma, utilizou-se de 366 pontos, aleatoriamente escolhidos no espaço, à semelhança do realizado por Gertner (2003). Os demais pontos, dos 626 totais, foram utilizados para comparar parâmetros de validação entre os métodos de krigagem ordinária (KO) e simulação sequencial gaussiana (SSG). O semivariograma foi estabelecido como a representação gráfica da semivariância estimada através da seguinte expressão: 29 ŷ(h) = N(h)2 1      N(h) 1i 2 ii h xZ)(x Z (31) em que, ŷ(h) é a semivariância estimada, e N (h) é o número de pares de valores observados Z (xi), Z (xi + h), separados pela distância h. O ajuste dos dados, a partir do semivariograma, possibilitou definir os seguintes parâmetros: efeito pepita (C0) – valor da semivariância quando h = 0; alcance da dependência espacial (a) – valor de h quando a semivariância se estabiliza próximo a um valor constante; patamar (C1 + C0) – valor da semivariância quando se obtém um valor constante próximo à variância dos dados. Foram ajustados aos dados os seguintes modelos: (a) esférico (Esf), ŷ (h) = C0 + C1 [1,5 (h/a) - 0,5 (h/a)3] para 0 < h < a e ŷ(h) = C0 + C1 para h > a; (b) exponencial (Exp), ŷ(h) = C0 + C1 [1 - exp (-3h/a)] para 0 < h < d, em que d é a máxima distância na qual o semivariograma é definido; (c) gaussiano (Gau), ŷ(h) = C0 + C1[1 - exp (-3h2/a2)]. A razão entre o efeito pepita e o patamar C0/(C0+C1), expressa em porcentagem, permitiu a classificação do grau de dependência espacial, a qual, segundo Cambardella et al. (1994), é considerada forte se esta razão for ≤ 25% do patamar, moderada quando ela está entre 25 e 75%, e fraca se a razão for > 75%. Para se determinar a existência ou não da dependência espacial, ajustaram-se semivariogramas, através do programa GS+ (GAMA DESIGN SOFTWARE, 2004). Para dirimir dúvida na escolha do melhor modelo, para um mesmo semivariograma, adotaram-se como critérios o melhor coeficiente de determinação (R2) e a menor soma de quadrados do resíduo (SQR). Para a elaboração dos mapas da distribuição espacial da erodibilidade em entressulcos utilizou-se de duas técnicas: 1) krigagem ordinária (KO), e 2) simulação sequencial gaussiana (SSG). Na krigagem ordinária, o modelo de estimação considerou as flutuações locais da média, limitando o domínio de sua estacionariedade para a vizinhança do local W(u). A média foi assumida ser constante, mas desconhecida. O estimador linear foi definido como (GERTNER, 2003): 30 Z*OK(u) =   )u(n 1 K(u) z(u) com   )u(n 1 OK(u) = 1 (32) em que, os n(u) pesos OK(u) são determinados de tal maneira que a variância do erro seja mínima. No caso da SSG, esta consistiu da aplicação de princípios da simulação sequencial para uma função aleatória com distribuição gaussiana. Considerou-se a área do estudo (A) como sendo composta pelos N nós de uma malha regular, tal que realizações {z(l)(u), u  A, com l=1,…,L}, representassem as L possibilidades da distribuição espacial dos valores do atributo z(u) ao longo da área A, em que cada realização foi denotada com o índice (l). Para tanto, uma função de distribuição acumulada condicionada (fdac) às variáveis aleatórias mais correlacionadas foi utilizada para obter valores z(u), para a variável aleatória z, em cada posição u  Α. Modelou-se a fdac considerando-se uma distribuição acumulada conjunta de m variáveis aleatórias, F(z1....zm (n)), condicionadas a um conjunto de n amostras do atributo: F(z1.... F(z1....zm (n)), = Prob{Zi  zi, i = 1 .... m (n)} (33) O processo descrito foi repetido nos N nós em diferentes sequências, ou seja, num caminho totalmente aleatório. Na prática, a fdac foi determinada pela média e variância obtidas por krigagem ordinária (KO), como descrito por King (2000). Na implementação da simulação sequencial condicionada, para obtenção de valores em m posições não-amostradas, seguiram-se as seguintes etapas: 1. Transformação dos z dados da distribuição das amostras disponíveis para uma distribuição gaussiana e a modelagem do variograma desses dados normalizados; 2. Definição de um caminho aleatório capaz de passar por todos os pontos; 3. Em cada ponto visitado, construir a função de distribuição condicional para o ponto, baseada nos pontos vizinhos com as amostras originais e previamente simuladas. Utilizar a KO para determinar a média e a variância da fdac gaussiana em cada nó; 4. Sorteio de um valor da fdac; 5. Adição do valor simulado ao conjunto de amostras, mais os pontos simulados anteriormente; 31 6. Repetição dos procedimentos de 1 a 4 para todos os pontos; 7.Transporte dos valores simulados dentro do espaço original, permitindo confrontar o conjunto de realizações com o modelo de histograma dos dados; 8. Repetição dos procedimentos de 1 a 7 para as L realizações. Para executar a SSG, como descrito anteriormente, utilizou-se o programa GS+ (GAMA DESIGN SOFTWARE, 2004) para um número de realizações L igual a 1.000. 32 4 RESULTADOS E DISCUSSÃO A quantidade média de sedimento erodido estimada com o modelo WEPP foi de 3097,63 kg ha-1. Esta quantidade é inferior à média anual estimada por Sanchez et al. (2009), com a Equação Universal de Perdas de Solo, para o Argissolo em Tabapuã (SP), que foi de 8260,00 kg ha-1. Martins Filho et al. (2009) obtiveram em parcelas experimentais, sob chuva simulada para o mesmo Argissolo com cultivo de cana-de-açúcar, perdas médias da ordem 2836,67 kg ha-1. Este último valor é próximo do estimado com o modelo WEPP. Os teores médios de fósforo disponíveis no solo e no sedimento erodido (Psolo e Psed) correspondem a classe “médio” (16-40 mg dm-3) para a cultura da cana-de- açúcar (VAN RAIJ et al., 1996). Tais teores são semelhantes aos encontrados em trabalhos realizados em solos sob cultivo da cana-de-açúcar (MARQUES JÚNIOR et al., 2008; CAMARGO et al., 2013a; OLIVEIRA et al., 2013). Em média, os sedimentos erodidos apresentaram-se 1,51 vezes mais enriquecidos por P que o solo original. Silva et al. (2012) obtiveram perdas médias de P disponível, em Latossolo sob cultivo de cana-de-açúcar, da ordem de 26 mg dm-3 para uma taxa de enriquecimento de 1,28. Tais resultados são inferiores aos obtidos no presente trabalho, mas são coerentes com a natureza mais erodível do Argissolo analisado. Martins Filho et al. (2009) obtiveram em parcelas experimentais (35 m2) taxas de enriquecimento do sedimento erodido por P de 2,7 a 3,8, as quais são superiores as avaliadas neste trabalho (Tabela 2). É possível observar na Tabela 2 que, os valores de taxa de enriquecimento WER obtida com o WEPP e com a função ajustada Ln(ER) = 2,682 – 0,278 Ln(Sed) (R2 = 0,69, p < 0,02) foram em média significativamente maiores que os obtidos com a equação Ln(ER) = 2 - 0,2 Ln(Sed), como apresentada por Perez-Bidegain et al. (2010). 33 Tabela 2. Valores da quantidade de sedimento erodido, concentrações de fósforo (P) no solo e no sedimento e taxas de enriquecimento por P para 17 encostas na área de estudo. Encosta Sed PSolo PSed ER PSed WER ER Pest (2) ER Pest (3) kg ha -1 mg dm -3 mg dm -3 H1 3978,68 63,9 81,5 1,28 1,51 1,47 1,41 H8* 4648,20 15,0 21,3 1,42 1,63 Ne 1,36 H10 3972,87 47,5 70,4 1,48 1,74 1,47 1,41 H15* 5001,56 39,9 55,6 1,40 1,56 Ne 1,35 H20 5329.75 32,5 46,9 1,44 1,65 1,36 1,33 H28 3972.40 10,7 17,1 1,60 1,56 1,47 1,41 H29* 2215,93 17,4 33,1 1,91 1,57 Ne 1,58 H33* 2214,47 37,1 69,9 1,88 1,72 Ne 1,58 H39* 2219,01 16,0 25,5 1,60 1,63 Ne 1,58 H44 2068,97 19,0 28,3 1,49 1,75 1,77 1,60 H45 1973,94 9,3 13,8 1,48 1,81 1,79 1,62 H52 2214,17 13,3 18,6 1,39 1,74 1,73 1,58 H54* 2213,40 17,3 27,5 1,58 1,56 Ne 1,58 H56* 3974,87 25,0 36,1 1,45 1,62 Ne 1,41 H70 2235,00 19,2 26,8 1,40 1,56 1,73 1,58 H88 2214,73 15,6 22,9 1,47 1,72 1,73 1,58 H89 2211,80 25,5 36,4 1,43 1,64 1,73 1,58 Médias 3097,63 24,6 37,2 1,51 1,65 a 1,63 a 1,50 b Sed é quantidade de sedimento erodido estimado com o WEPP derivado de cada encosta (kg ha -1 ), PSolo é o conteúdo de fósforo no solo (kg kg -1 ), PSed é quantidade de P no sedimento (kg ha -1 ), ER é a taxa de enriquecimento por fósforo, (2) é a ER obtida com a função ajustada Ln(ER) = 2,682 – 0,278 ln(Sed) (R 2 = 0,69, p < 0,02) obtida com os dados das encostas grafadas com *, Ne representa valor não estimado, (3)é a ER obtida com a função Ln(ER) = 2 - 0,2 Ln(Sed). Médias seguidas de mesma letra na linha não diferem entre si pelo teste t (p < 0,05). Os valores do coeficiente Nash-Sutcliffe indicam o ajustamento dos dados simulados aos observados em relação à reta 1:1, podendo variar de -∞ a 1. O NS está associado à eficiência de estimativa de enriquecimento por fósforo disponível do sedimento erodido. Como descrito por Viola et al. (2012), tem-se que: NS > 0,65 o modelo é considerado muito bom; 0,54 < NS < 0,65 o modelo é considerado bom e entre 0,5 e 0,54, satisfatório. Nos casos da Tabela 3 os três modelos testados foram classificados como muito bons. Na análise do erro quadrático médio (RMSE) observa-se que o valor com o modelo 1 foi o menor em relação aos modelos 2 e 3 (Tabela 3). Valores de RMSE iguais a 0 são considerados como iguais. 34 Os valores de MD > 0 indicam que os valores de P estimados pelo modelo superestimam os valores observados ou medidos, o que é o caso dos modelos 2 e 3. Já o modelo 1 tende a subestimar os valores observados com as estimativas. Quanto ao índice de concordância (d) de Wilmott, (1982), o qual expressa exatidão entre valores estimados e observados, verificou-se que os três modelos são satisfatórios (d > 0,75), segundo critérios adotados por Lopes et al. (2014). Pode-se afirmar que os modelos proporcionaram alto grau de exatidão, uma vez que os seus índices de concordância (d) foram próximos da unidade. Tabela 3. Desempenho de três métodos de estimativa da taxa de enriquecimento por fósforo do sedimento erodido. Modelo NS RMSE MD d 1. WER 0,937 0,024 0,020 0,995 2. ln(ER) = 2,682 – 0,278 ln(Sed) 0,969 0,019 0,010 0,997 3. ln(ER) = 2 - 0,2 ln(Sed) 0,981 0,013 -0,002 0,998 ER é a taxa de enriquecimento por fósforo estimada pelos modelos 1 e/ou 2, WER é a taxa de enriquecimento por fósforo estimada com o modelo WEPP - modelo 1, NS é o coeficiente de Nash-Sutcliffe, RMSE é o erro quadrático médio, MD é a diferença média e d o índice de concordância. Na Figura 4, os valores estimados em função dos observados de fósforo disponível no sedimento erodido são apresentados em relação à reta 1:1. Na Figura 4 ficam evidentes as tendências dos valores estimados pelos modelos 1 e 2 superestimarem os valores observados de P no sedimento erodido. O contrário pode ser constatado com o modelo 3. Tais fatos foram confirmados pelos valores positivo e negativo de MD (Tabela 3). 35 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Fósforo observado, Pobs, kg ha-1 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 F ó s fo ro e s ti m a d o , P e s t, k g h a -1 1:1 Pest (1) = -0,007 + 1,167 Pobs Pest (2) = -0,0196 + 1,061 Pobs Pest (3) = -0,0007 + 0,991 Pobs Figura 4. Valores medidos e estimados para fósforo disponível no sedimento erodido e de taxas de enriquecimento (WER e ER) obtidos como: (1) WER do modelo WEPP; (2) Ln(ER) = 2,682 – 0,278 Ln(Sed); (3) Ln(ER) = 2 - 0,2 Ln(Sed). Os resultados apresentados corroboram com os trabalhos realizados por Perez-Bidegan et al. (2010) e Elliot et al. (2015), no sentido de ser possível incorporar ao WEPP um algoritmo ou rotina para predizer o fósforo derivado de encostas com o sedimento erodido. Deste modo, foi estimado para a área da bacia composta pelas 89 encostas a distribuição espacial das perdas de P disponível com o sedimento erodido. Para tanto, os modelos apresentados na Tabela 3 para estimar a taxa de enriquecimento por fósforo do sedimento erodido foram utilizados. É possível verificar que as taxas de enriquecimento WER não diferiram significativamente entre encostas convexas e côncavas (Tabela 4). Contudo, as taxas de enriquecimento nas encostas côncavas foram significativamente maiores que as das convexas, quando obtidas com os métodos (2) e (3). Quando se analisou o P disponível originalmente no solo, não se constatou diferença significativa entre as formas convexas e côncavas. Já no sedimento e na enxurrada verificou-se que as concentrações médias de P disponível foram 36 significativamente maiores nas encostas côncavas em relação às convexas (Tabela 4). Tabela 4. Fósforo no solo, Taxas de enriquecimento, fósforo no sedimento e na água da enxurrada, e variáveis da erosão do solo em função da curvatura do relevo. Análises Variáveis Curvatura do relevo Solo PSolo, mg dm -3 Convexa Côncava 32,61 a 36,02 a WER (1) 1,62 a 1,63 a Enriquecimento ER (2) 1,57 a 1,39 b ER (3) 1,70 a 1,44 b P (1), kg ha -1 0,13 b 0,24 a Sedimento P (2), kg ha -1 0,13 b 0,22 a P (3), kg ha -1 0,14 b 0,22 a P (1), mg L -1 0,07 b 0,11 a Enxurrada P (2), mg L -1 0,07 b 0,10 a P (3), mg L -1 0,07 b 0,10 a A, kg ha -1 ano -1 2397,45 b 4344,67 a Erosão CS, kg L -1 0,00132 b 0,00197 a VE, m 3 ano -1 3144,37 b 6312,48 a Produção de Sedimento, kg ano -1 4617,67 b 14153,54 a Deposição de sedimento, kg ano -1 161,49 b 1094,20 a Taxas de enriquecimento (WER e ER); P(1), P(2) e P(3) são as concentrações de P disponível no sedimento e na enxurrada, para valores de taxas de enriquecimento (WER e ER) obtidos com: (1) WER do modelo WEPP; (2) Ln(ER) = 2,682 – 0,278 Ln(Sed); (3) Ln(ER) = 2 - 0,2 Ln(Sed), respectivamente; A é a taxa de erosão; CS é concentração de sedimento na enxurrada; VE é o volume de enxurrada. Médias seguidas de mesma letra na linha não diferem entre si pelo teste t (p < 0,05). Todas as médias das variáveis relativas à erosão do solo, tais como taxa de erosão líquida (A), concentração de sedimentos na enxurrada (CS), volume de enxurrada (VE), produção e deposição de sedimentos foram significativamente maiores nas encostas côncavas em comparação com as convexas (Tabela 4). Estes resultados contradizem os de Sanchez et al. (2009), os quais foram obtidos na mesma área do presente estudo. Sanchez et al. (2009) estimaram, com a Equação Universal de Perdas de Solo (USLE), maiores perdas de solo por erosão para as pedoformas convexas em relação às côncavas. Numa outra pesquisa, nesta mesma área de estudo, Ibiapina (2015) demonstrou que o modelo WEPP propiciou maior precisão e exatidão na estimativa 37 dos valores de erosão em relação à USLE. Os resultados expressos na Tabela 4 são justificáveis por resultados obtidos por Camargo et al. (2013a) e Barbosa (2015). Camargo et al. (2013a) avaliaram a mineralogia do Argissolo utilizado no presente trabalho, em curvaturas do relevo, os quais concluíram que os maiores teores de hematita e goethita foram encontrados na área convexa e que há predominância da hematita em ambas curvaturas (convexa e côncava). Na área convexa a hematita apresentou o menor diâmetro médio do cristal. Foi verificado haver predominância de caulinita na área côncava. A presença de caulinita implica no desenvolvimento de macroestrutura do tipo em blocos, como observado por Ferreira et al. (1999), devendo originar solo com maior densidade, maior proporção de poros pequenos e menor permeabilidade. Logo, a erodibilidade deverá ser maior em solos cauliníticos. A expressão de tais fatos são os maiores valores de A (Figura 5), CS, VE, produção e deposição de sedimentos significativamente maiores nas encostas côncavas em relação às convexas. Figura 5. Padrão espacial obtido por krigagem ordinária da erosão líquida estimada com o modelo WEPP. 38 Para o Argissolo Vermelho-Amarelo analisados neste trabalho, Barbosa (2015) demonstrou que os valores de erodibilidade em sulcos (Kr) nas posições de forma convexa são inferiores aos de forma côncava, 0,0053 kg N-1 s-1 e 0,0146 kg N- 1 s-1, respectivamente. Já os valores de tensão cisalhante crítica são de 2,1887 N m- 2 na forma convexa e de 2,7260 N m-2 na forma côncava. Tais resultados são coerentes com os padrões de erosão estabelecidos pelo modelo WEPP na bacia estudada, obtidos por Martins Filho (2007) e, ainda, corroboram os resultados de Camargo et al. (2013a). A mineralogia local tem influência nos atributos do solo e na hidrologia local, o que afeta o grau de erosão local. Assim, tal fato também influenciará a adsorção de P ao sedimento erodido. Em função dos resultados, é possível averiguar que as maiores perdas de P ocorrem em encostas de formas côncavas que nas convexas. A área convexa, segundo Camargo et al. (2013a), apresenta os maiores valores de diâmetro médio ponderado, diâmetro médio geométrico, agregados >2 mm, agregados 2-1 mm, volume total de poros, e teor de água no solo e os menores valores de densidade do solo e resistência à penetração nas duas profundidades estudadas em relação à área côncava. Camargo et al. (2013a) verificaram, ainda, que a área côncava apresentou um expressivo número de atributos do solo com maior variabilidade espacial. Um aspecto relevante é que a mineralogia da fração argila influencia na adsorção de P, como ressaltaram Camargo et al. (2013a). Fenômenos, como a adsorção de elementos no solo estão fortemente relacionados com a cristalinidade mineral. Assim, o P tende a ser adsorvido por goethita (Gt), com baixa cristalinidade e uma elevada SSA, bem como pela gibbsita (Gb). Isto permite a utilização da técnica descrita neste estudo em alternativas práticas de gestão, com vista a minimizar os efeitos da adsorção de fósforo no solo, por exemplo. Os valores de WER variaram de 0,37 a 2,27, com valor médio de 1,63. ER (1) e ER(2) variaram de 1,33 a 1,62 e de 1,35 a 1,77, nesta ordem (Tabela 5). Média e mediana não foram muito discrepantes para WER, ER(1) e ER(2). Valores estimados de P disponível variaram: P(1) de 0,006 a 4,050 kg ha-1; P(2) de 0,006 a 4,757 kg ha-1; P(3) de 0,007 a 4,929 kg ha-1. Já o as concentrações de P disponível na enxurrada variaram: P(1) de 0,004 a 1,751 mg L-1; P(2) de 0,003 a 2,057 mg L-1; 39 P(3) de 0,003 a 2,131 kg ha-1. O volume de enxurrada variou de 668,30 a 11084,70 m3 ano-1, com média de 4612,03 m3 ano-1. A erosão do solo 1973,74 a 5329,75 kg ha-1 ano-1, com média de 3299,52 kg ha-1 ano-1. Valores de concentração de sedimentos variaram de 0,001 a 0,003 kg L-1, com média de 0,002 kg L-1 (Tabela 5). Adotando-se a classificação de Warrick e Nielsen (1980), todos os valores de estimativas das concentrações de P disponível no sedimento e enxurrada, volume de enxurrada, erosão líquida e concentração de sedimento na enxurrada apresentaram valores altos de CV (CV ≥ 24%), com exceção dos valores de taxas de enriquecimento, que apresentaram baixo valor de CV (≤ 12%). Estes resultados refletem a variação textural, química e mineralógica da área de estudo, o que já foi avaliado por Martins Filho (2007), Sanchez et al. (2009) e Camargo et al. (2013a,b). É comum resultados relacionados à erosão do solo apresentarem altos coeficientes de variação e a não normalidade dos dados amostrais (Tabela 5). Tabela 5. Estatística descritiva das taxas de enriquecimento, fósforo no sedimento e na água da enxurrada, e variáveis da erosão do solo para 626 pontos nas encostas convexas e côncavas. Estatística WER ER (2) ER (3) P (1)Sed P (2)Sed P (3)Sed P (1)Enx P (2)Enx P (3)Enx VE A CS kg ha-1 kg ha-1 kg ha-1 kg L-1 kg L-1 kg L-1 m3 ano-1 kg ha-1 ano-1 kg L-1 Média 1,63 1,49 1,56 0,180 0,169 0,176 0,087 0,082 0,086 4612,03 3299,52 0,002 Mediana 1,61 1,58 1,71 0,071 0,064 0,069 0,034 0,031 0,033 3861,70 2225,34 0,001 Mínimo 0,37 1,33 1,35 0,007 0,006 0,007 0,004 0,003 0,003 668,30 1973,94 0,001 Máximo 2,27 1,62 1,77 4,050 4,757 4,929 1,751 2,057 2,131 11084,70 5329,75 0,003 D. Padrão 0,18 0,11 0,161 0,325 0,335 0,348 0,148 0,154 0,160 2479,21 1232,21 0,045 ₓ 10-2 Assim. 0,105 -0,211 -0,205 5,277 6,549 6,506 4,753 5,849 5,788 0,789 0,313 0,680 (0,32) (0,40) (0,36) (0,46) (0,41) (0,47) Kurtose 4,506 -1,839 -1,845 42,385 66,004 65,333 34,546 52,822 51,666 0,151 -1,700 -0,914 (-0,48) (-0,32) (-0,35) (-0,31) (-0,42) (-0,29) N 626 626 626 626 626 626 626 626 626 626 626 626 CV (%) 11,268 7,387 10,219 180,882 199,058 197,620 170,594 187,992 186,879 53,755 37,345 27,880 Taxas de enriquecimento (WER e ER); P(1), P(2) e P(3) são as concentrações de P disponível no sedimento (Sed) e na enxurrada (Enx), para valores de taxas de enriquecimento (WER e ER) obtidos com: (1) WER do modelo WEPP; (2) Ln(ER) = 2,682 – 0,278 Ln(Sed); (3) Ln(ER) = 2 - 0,2 Ln(Sed), respectivamente; VE é o volume de enxurrada; A é a taxa de erosão; CS é concentração de sedimento na enxurrada. ( ) = valores após transformação Log. A caracterização da estrutura da variabilidade espacial foi determinada com o ajuste de modelos teóricos aos variogramas experimentais (Tabela 6). Utilizou-se a transformação logarítmica nos dados de PSed e PEnx, visando à correção da assimetria da distribuição (Tabela 5). A transformação logarítmica é recomendada 40 para conjunto de valores com assimetria positiva maior que a unidade (OLIVER e WEBSTER, 2014). O modelo esférico apresentou os melhores ajustes para os atributos PSed, e PEnx e A, (Tabela 6). Bahia et al. (2015) também utilizaram o modelo esférico para atributos analisados e, somente para Hm, estes autores utilizaram o modelo exponencial, diferente deste trabalho. O modelo esférico se ajusta a atributos que apresentam variações abruptas ao longo da paisagem, e em ciência do solo é um dos mais utilizados (VIEIRA, 2000). Estas variações podem estar relacionadas aos tipos de material de origem (RAUCH, 2011), relevo (SIQUEIRA et al., 2010; CAMARGO et al., 2013a), solo (MONTANARI et al., 2012). Todos os atributos relativos às perdas de fósforo apresentaram moderado grau de dependência espacial (25% < GDE < 75%, Tabela 6), conforme classificação de CAMBARDELLA et al. (1994). Já para o atributo erosão do solo (A) obteve-se alto grau de dependência espacial (GDE < 25%). Os valores do GDE indicam que estes atributos têm maior influência dos compartimentos da paisagem, apresentando maior potencial para refinamento dos limites de precisão (BAHIA, 2016). Verificou-se que os atributos estudados apresentaram alcances de 424 m a 507 m para PSed e PEnx, enquanto que A apresentou 829 m de alcance. Alcances maiores indicam uma maior continuidade da distribuição espacial dos atributos. Oliveira et al. (2013) estabeleceram alcances para 349 m para P disponível no Argissolo da bacia em apreço. Tal diferença entre os alcances pode ser justificada, pois Oliveira et al. (2013) verificaram que áreas com altos teores de P disponível no solo exibiram distintos padrões espaciais, evidenciando menor continuidade espacial. Outro aspecto é que os alcances para P disponível, no presente trabalho, foram obtidos para fósforo no sedimento e na enxurrada, os quais estão estritamente ligados à erosão (A), a qual apresentou maior continuidade espacial demonstrada pelo alcance de 829 m. Os atributos apresentaram elevados coeficientes de determinação (R2), indicando que mais de 99% da variabilidade espacial são explicados pelos modelos ajustados (Tabela 6). 41 Tabela 6. Modelos e parâmetros dos semivariogramas simples ajustados aos dados. Transf. Atributo Modelo C0 C0+C1 GDE (%) a (m) R2 SQR Log P(1)Sed Esférico 0,865 1,731 50,0 507 0,99 7,679 ₓ 10-3 Log P(2)Sed Esférico 0,841 1,703 49,4 480 0,99 6,428 ₓ 10-3 Log P(3)Sed Esférico 0,852 1,732 49,2 478 0,99 6,572 ₓ 10-3 Log P(1)Enx Esférico 0,770 1,541 50,0 454 0,99 6,489 ₓ 10-3 Log P(2)Enx Esférico 0,774 1,549 50,0 433 0,99 9,237 ₓ 10-3 Log P(3)Enx Esférico 0,0778 1,557 50,0 424 0,99 9,237 ₓ 10-3 Normal A Esférico 120000 1955000 0,6 829 0,99 1,090 ₓ 10-2 Transf. = transformação; P(1), P(2) e P(3) são as concentrações de P disponível no sedimento (Sed) e na enxurrada (Enx), para valores de taxas de enriquecimento (WER e ER) obtidos com: (1) WER do modelo WEPP; (2) Ln(ER) = 2,682 – 0,278 Ln(Sed); (3) Ln(ER) = 2 - 0,2 Ln(Sed), respectivamente; A é a taxa de erosão (kg ha -1 ano -1 ). Na Figura 6 pode-se observar que não há grande diferença no padrão espacial das perdas de P disponível na área da Bacia, com os três modelos testados para estimar a taxa de enriquecimento (ER ou WER). No caso da Figura 6a a taxa de enriquecimento foi a estimada pelo WEPP, ou seja, a WER. Já no caso das Figuras 6b e 6c, a ER foi estimada com modelos empíricos. As concentrações médias anuais de P disponível ou solúvel na enxurrada, derivadas em 17 encostas da Bacia, foram da ordem de 0,059 a 0,067 mg L-1 (Tabela 7). Tais valores são preocupantes, pois são três vezes maiores que o limite crítico de fósforo total presente na água que é 0,020 mg L-1. Logo, torna-se relevante utilizar modelos que não estimem apenas o escoamento superficial e o sedimento derivado das encostas, mas que também avaliem as perdas de P com o sedimento e água da enxurrada, como mencionado e realizado por Elliot et al. (2015). 42 a Kg ha -1 b c Figura 6. Distribuição espacial obtida por simulação sequencial gaussiana das perdas de fósforo disponível no sedimento erodido, com valores de taxas de enriquecimento (WER e ER) obtidos como: (a) WER do modelo WEPP; (b) Ln(ER) = 2,682 – 0,278 Ln(Sed); (c) Ln(ER) = 2 - 0,2 Ln(Sed). 43 Tabela 7. Concentrações do sedimento erodido, fósforo no solo e na água da enxurrada em 17 encostas da bacia. Encosta CS PSolo P (1) P (2) P (3) kg/L mg dm-3 mg L-1 mg L-1 mg L-1 H1 0.001718 63.9 0.166 0.154 0.160 H8 0.002076 15.0 0.051 0.043 0.044 H10 0.001718 47.5 0.142 0.115 0.119 H15 0.002261 39.9 0.140 0.121 0.123 H20 0.002657 32.5 0.143 0.115 0.116 H28 0.001718 10.7 0.029 0.026 0.027 H29 0.001234 17.4 0.034 0.034 0.037 H33 0.001234 37.1 0.079 0.073 0.079 H39 0.001234 16.0 0.032 0.031 0.034 H44 0.001256 19.0 0.042 0.038 0.042 H45 0.001271 9.3 0.022 0.019 0.021 H52 0.001234 13.3 0.029 0.026 0.028 H54 0.001234 17.3 0.033 0.034 0.037 H56 0.001718 25.0 0.070 0.060 0.063 H70 0.001234 19.2 0.037 0.037 0.041 H88 0.001234 15.6 0.033 0.030 0.033 H89 0.001234 25.5 0.052 0.050 0.054 Média 0.001545 24.9 0.067a 0.059a 0.062ª CS é concentração de sedimento na enxurrada; Psolo é a concentração de fósforo disponível no solo; P(1), P(2) e P(3) são as concentrações de P na enxurrada, com valores de taxas de enriquecimento (WER e ER) obtidos com: (1) WER do modelo WEPP; (2) Ln(ER) = 2,682 – 0,278 Ln(Sed); (3) Ln(ER) = 2 - 0,2 Ln(Sed), respectivamente. A distribuição espacial da probabilidade de ocorrer uma concentração de P > 0,020 mg L-1 na enxurrada, considerando-se os três métodos de determinação da taxa de enriquecimento, pode ser observada na Figura 7. Verificou-se que o padrão de distribuição espacial de concentrações de P > 0,020 mg L-1 são semelhantes, entre os modelo1, 2 e 3 (Figuras 7a, 7b e 7c). Para estes modelos, nas regiões mais críticas (p>0,75), a probabilidade de P > 0,020 mg L-1 ocorre em mais de 81% da área da bacia (Tabela 8). Para probabilidades (p) <0,25, 0,25 a 0,50 e de 0,50 a 0,75, as percentagem de áreas representativas ocupadas por estas classes na bacia são em média 1,7%, 6,6% e 10%, respectivamente. 44 a Probabilidade p (P > 0,02 mg/L) b c Figura 7. Distribuição espacial obtida por simulação sequencial gaussiana da probabilidade da concentração de fósforo disponível ser maior que 0,020 mg L-1 na enxurrada, com valores de taxas de enriquecimento (WER e ER) obtidos como: (a) WER do modelo WEPP; (b) Ln(ER) = 2,682 – 0,278 Ln(Sed); (c) Ln(ER) = 2 - 0,2 Ln(Sed). 45 Tabela 8. Percentagem da área da bacia por intervalo de probabilidade (p) para a concentração de fósforo disponível ser maior que 0,020 mg L-1 na enxurrada. Intervalo probabilidade (p) Modelo1 (a) Modelo 2 (b) Modelo 3 (c) ---------------------------- % da Área da bacia ----------------------------- 0 a 0,25 1,6 1,9 1,6 0,25 a 0,50 6,1 7,5 6,1 0,50 a 0,75 10,3 9,7 9,9 >0,75 82,1 80,9 82,3 Valores