A utilização de modelos matemáticos na análise de sistemas de recursos hídricos é bastante discutida na literatura, entretanto todos apresentam simplificações quanto a representação matemática da dinâmica dos seus processos ou nas aproximações numéricas das equações não-lineares que descrevem o seu comportamento. Dependendo da aplicabilidade, tais simplificações podem comprometer uma tomada de decisão. Em geral, os modelos de simulação disponíveis para planejamento de sistemas de recursos hídricos a nível de bacia hidrográfica, apesar de serem bastante versáteis, não conseguem considerar certas complexidades devido a algumas limitações matemáticas dos seus algoritmos. Neste contexto foi proposto um modelo de simulação para esta finalidade, que tem como principais características o desenvolvimento de um modelo quali-quantitativo integrado com os múltiplos usos de águas superficiais, a incorporação de não-linearidades dos processos hidráulicos e operacionais e uma função multiobjetivo que otimize, segundo critérios de prioridade, o atendimento das demandas quantitativas hídricas dos diferentes usos e qualitativas relativas ao atendimento de metas de enquadramento para uso, segundo as normas do CONAMA, assim como o atendimento de outras metas operacionais. Os parâmetros de qualidade da água envolvidos são: Demanda Bioquímica de Oxigênio, Oxigênio Dissolvido, Fósforo Total, Nitrogênio Total, Clorofila-a e Coliformes Fecais. As limitações operacionais e hidráulicas foram incluídas no modelo através de restrições lineares (balanço hídrico nos nós, volumes metas, volumes mínimos, capacidade dos componentes, etc.) e não-lineares (vazão vertida máxima, vazão descarregada máxima, área da superfície líquida do reservatório, etc.). As restrições não-lineares foram tratadas a partir de aproximações por segmentos lineares (artifícios de linearização) e resolvidas num processo iterativo (programação linear sequencial) até a convergência do processo a uma tolerância desejada para o erro relativo da função objetivo. As equações do balanço de massa para os parâmetros de qualidade da água, que geralmente requerem o produto de variáveis de decisão, foram linearizadas pelo Método das Aproximações Lineares e agregadas, num mesmo algoritmo, às equações do balanço hídrico. Para avaliar o desempenho do atendimento das demandas hídricas, foram incluídos, no modelo, alguns indicadores de análise de desempenho, como a confiabilidade, vulnerabilidade, resiliência e sustentabilidade. Sua aplicabilidade foi verificada através de uma análise sistêmica dos reservatórios Engenheiro Ávidos e São Gonçalo da bacia hidrográfica do Rio Alto Piranhas – PB e suas demandas quali-quantitativas. A análise deste sistema foi realizada através da idealização de quatro cenários de operação que envolvia o uso de prioridades diferentes e a inclusão de uma Estação de Tratamento de Esgoto – ETE no último cenário. O horizonte de tempo foi de 360 meses. Foram analisados o atendimento às demandas e o enquadramento dos níveis de concentração dos parâmetros de qualidade da água à Classe II do CONAMA 357/05. Os resultados demonstraram que todas as restrições do sistema foram satisfeitas e que as demandas quantitativas obtiveram diferentes garantias de atendimento de acordo com o cenário utilizado. A inclusão da ETE considerada no cenário 04 foi imprescindível para aumentar o atendimento aos usos da água e melhorar a qualidade da água no sistema. O modelo proposto mostrou-se eficiente, quer seja quanto ao tempo de processamento quanto ao atendimento de todas as restrições impostas ao sistema. Apresentou, também, um erro numérico relativamente baixo.