This paper systematically simulates how permeability is influenced by open vugs of different sizes, shapes, and abundance in three-dimensional reservoir volumes. A total of 1920 simulations analyzed using logistic regression reveal mathematical relationships that predict the probability (P) of high permeabilities from connected passageways of vugs, given known vug shape (Ar), size (Si), and abundance (Ab): P = 1/Exp(−24.139 + [0.702 × Ab] + [2.857 × Ar] + [0.425 × Si]) for oblate vug shapes; and P = 1/Exp(−28.803 + [0.973 × Ab] + [3.147 × Ar] + [0.081 × Si]) for prolate vug shapes. Standardized regression coefficients of simulation results show that vug abundance has the greatest impact on vug connectivity, followed by vug shape, and then vug size. Anisotropy of permeability can be predicted from the volumetric abundance of the largest network of touching vugs: omnidirectional permeability at >5%; transverse permeability at ≥ 2% but ≤ 5%; transverse permeability or isolated vugs without connections across the volume at 0.4% to 2%; and isolated vugs without connections across the volume below 0.4%. The workflow for calculating the probability that vugs will connect to form permeability passageways, represented as a log, are (1) use well logs, such as nuclear magnetic resonance, to calculate vug size and abundance; (2) determine vug shape by other means; and (3) input data into logistic regression algorithms. This workflow was applied using a well-studied subsurface example from the Arbuckle Group in Kansas. The application verifies that this new workflow (and logistic regression equations) provides a successful approach for identifying vuggy permeability zones and their anisotropies in the subsurface (A PCT [Patent Cooperation Treaty] International Patent Application publication number WO2022132784 titled [Algorithms for Predicting Vug Connectivity] was filed on December 14, 2021, related to this work).