First it would be good to simplify the problem - you should test the method first with some values before potting it!
julia> severinghaus(1u"Torr")
ERROR: DimensionError: 1 Torr^3 and 150 Torr are not dimensionally compatible
In this case the severinghaus function is kind of an empirical hack without units. Pressure cubed plus pressure, which isn’t a thing Unitful.jl can do.
But these functions exist and can be useful, so I usually do the opposite to what you have done and strip the units out in this situation. Stripping with u"Torr" makes sure its the right units fed into the equation, instead of some other Pressure type:
severinghaus(ppO2) = (23400 * (ppO2^3 + 150 * ppO2)^-1 + 1)^-1
severinghaus(ppO2::Unitful.Pressure) = severinghaus(Unitful.ustrip(u"Torr", ppO2))
plot(severinghaus, (1:1:150)u"Torr")
Also notice you need to pass a second argument to plot with unitful values - these are what is sent to the servinghouse function, plots wont actually know to convert that from just passing unitful xlims.
If the output actually has a unit, you can add it back at the end of the unitful aware severinghaus method
Edit: after discussions below, I think you should put units on all the constants instead, even though they probably aren’t in the original.