How do I set the proper boundary conditions in solve_bvp?

bc : callable
     Function evaluating residuals of the boundary conditions [← emphasis mine].
     The calling signature is bc(ya, yb), or bc(ya, yb, p) if parameters are present.
     All arguments are ndarray: ya and yb with shape (n,), and p with shape (k,).
     The return value must be an array with shape (n + k,).

The return value is the violation of the boundary condition, or in other words you need a function that evaluates to zero when the boundary conditions ares satisfied.

def bc(ya, yb): return ya[0]-23.8, yb[1]

Well, I hesitated before posting this answer as I would have rather checked the solution more thoroughly but as I see that this is still open here all what I have found. See the code below but some explanations first:

• first I changed your ODEs so its consistent with your ODE in your question
• I have added a modified right-hand side (fun) adding parameters - this is to cover the now changed boundary conditions from transforming the ODE. I kept close to the scipy tutorial so the explanations match up. You might want to try different starting values for the parameter p.
• the same for the boundary conditions (bc1) - the logic in naming the bc parameters is Sn: n=0 boundary a, n=1 boundary b; [0] = y, [1] = y' - so the name covers the side of the boundary, the index the degree of derivative
• since I do not know the physics of the problem I added a graph to show the solution. This seems to look right but mind the infinity-condition: this is now set to 0 at the b-boundary (the "other side"). Infinity can not be covered, the only solution I see is to expand that limit outward.

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

# constants
D = 200.6
v = 0.02
k = 0.20

# split second order ODE into 2 first order ODEs
'''def sulfate_profile(U, z):
    print(U, z)
    dUdz = [U[1], (v*U[1] + U[0]/(U[0]+k))/D]
    print(dUdz)
    return dUdz'''

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html
def sulfate_profile(z, S): # compare to doc: x->z, y->S
    return np.vstack((S[1], v/D*S[1] + S[0]/(S[0]+k)/D))

def fun(z, S, p):
    k0, k1 = p[0], p[1]
    return np.vstack((S[1], k0*S[1] + k1*S[0]/(S[0]+k)))

# I initiated these boundary conditions but I don't really know what they mean...
def bc(S0, S1):
    return np.array([S0[0], S1[0]])

def bc1(S0, S1, p):
    k0 = p[0]
    return np.array([S0[0], S1[0], S0[1] - k0*23.8/(23.8+k), S1[0]]) # Sn: n=0 boundary a, n=1 boundary b; [0] = y, [1] = y',...

z_list = np.linspace(0, 100, 100)
S_a = np.zeros((2, z_list.size))
S_a[0, 0] = 23.8
S_a[0, 99] = 0 # b-boundary which is inf - expand here

# solve SO4 concentration ODE over depth range (z_list)
#sulfate_reduction = solve_bvp(sulfate_profile, bc, z_list, S_a)
sulfate_reduction = solve_bvp(fun, bc1, z_list, S_a, p=[1,1])

#print(sulfate_reduction)
x_plot = np.linspace(0, 1, 100)
y_plot_a = sulfate_reduction.sol(x_plot)[0]
plt.plot(x_plot, y_plot_a, label='sulfate_b_dpth')
plt.show()

This results in

which looks like a depth distribution (if that is what it is).