We developed a methodology to estimate maximum brine injection rates in subsurface formations across wide geographic areas using inverse modeling-based optimization techniques. We first defined geographic areas where groundwater was too saline to meet the standard for drinking water and where sufficient confining units existed above and below the injection layers. We then assumed concurrent brine injection into a system of wells on a consistent 25 km x 25 km spacing across the entire modeled area. Taking advantage of symmetry, we represented each 25 km x 25 km injection area as a 12.5 km-long one-dimensional radial model, divided into 100 logarithmically-sized grid blocks. A single layer of grid blocks was used because homogenous porous media were assumed. Brine injection was simulated into the leftmost (innner) grid block, and the injection rate was automatically adjusted to meet a maximum pressure buildup to 80% of the fracturing pressure, estimated as the least principal stress, at the injection location. A secondary constraint of 1 bar maximum pressure increase at the right-most (far-field boundary) grid block after 50 years of injection was applied. We demonstrated this method on three stratigraphic layers that overlie the Mt. Simon Sandstone (MSS) in the Illinois Basin, as well as in the MSS itself, because the MSS is a well-known CO2 injection target with a large estimated CO2 storage capacity. CO2 storage in the MSS could be optimized by extracting brine from that formation and injecting it elsewhere, so the brine injection rates estimated with the models contained herein could help to refine CO2 storage capacity estimates.