using NonlinearSolve
using ModelingToolkit
# !!! used UNITS !!!
# eq/liter
# Pa
# constants
mmHG = 133.322 # Pa
Kc = 2.46e-11/mmHG
Ka = 3e-7
@variables HA A H HCO3
@parameters pCO2 Atot SID
eqs = [
H * A ~ Ka * HA
H * HCO3 ~ Kc * pCO2
Atot ~ HA + A
SID ~ HCO3 + A
]
@mtkbuild ns_stewart = NonlinearSystem(eqs)
display(ns_stewart)
println()
# Initial values (guesses)
guesses_stewart = (
HCO3 => 24.9e-3,
H => 3.95e-8,
A => 17.1e-3,
HA => 2e-3
)
ps = (
Atot => 19e-3,
SID => 41.7e-3,
pCO2 => 5.3e3
)
prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)
sol = solve(prob)
@show sol
@show H = -log10(sol[1])
println("HCO3 = $(sol[2]*1000) meq/liter")

NEW Results:

julia> include("stewart1.jl")
Model ns_stewart with 2 equations
Unknowns (2):
H
HCO3
Parameters (3):
pCO2
Atot
SID
sol = [3.9275240698625016e-8, 0.02489948137605699]
H = -(log10(sol[1])) = 7.405881144702646
HCO3 = 24.89948137605699 meq/liter

My interpretation:
In equ1 I used the equations I found in the article you referred. The variable H here is a part of the LHS expression.
In equ2 the variable H is explicit mentioned in 2 equations. For that reason H was eliminated by MTK.

In one case H is nonlinearly related to another variable. It’s just a case the tearing algorithm does not handle yet. In fact, it’s hard to reduce because it’s not generally valid to divide like that because in theory that unknown could be zero.