Вулканические извержения и сопровождающие их потоки лавы представляют собой значительную опасность для населения, построек и инфраструктуры региона. Лава может занимать большие пространственные области, для которых детальное трехмерное моделирование процесса ее течения сводится к решению дискретных задач очень большой размерности и не всегда оказывается эффективным. При достаточно малом отношении вертикального размера потока к его горизонтальному размеру применяются математические модели, основанные на усредненных по глубине уравнениях движения вязкой среды. В данном исследовании такая модель состоит из уравнений для глубины лавы, двумерных уравнений ее движения, кинетики кристаллов, уравнения теплового баланса, которое учитывает нелинейный конвективный и радиационный обмен энергии с внешней средой, энергию диссипации и скрытую теплоту кристаллизации. Математическая модель реализована численно в пакете OpenFOAM с открытым исходным кодом. Пакет позволяет использовать для осуществления вычислительных экспериментов современные высокопроизводительные кластеры и адаптировать задачу к конкретным физическим аспектам моделируемого природного процесса. Проведена верификация кодов путем сравнения аналитического решения задачи с решением согласно модели с уравнениями, описывающими движение в пространственной области двухфазной несжимаемой жидкости. Исследовано влияние на поток реологических характеристик на примере представления его моделью Ньютона по сравнению с моделью Бингама и нелинейной моделью Гершеля–Балкли. Нелинейная реология рассматриваемой жидкости учитывает зависимости фактической вязкости лавового потока от температуры, скорости сдвига, предела текучести (при этом предел текучести и степень нелинейности являются функциями температуры. Параллельные компьютерные коды реализованы с помощью интерфейса OpenMPI на вычислительных кластерах с общей и распределенной памятью под управлением ОС Linux. Проведено профилирование кодов для многоядерных процессоров с общей памятью.