Just two disclaimers to start with:
I'm not using it presently.
Due to this sentence and the lack of broader context in your question (e.g., if there are any other stochastic variables in your model besides $d_i$), I assume that right now you have a completely deterministic model and this is the first stochastic description you want to add to the model.
Suppose I know that di is the center of some normal distribution, and I also know vi, the variance for that distribution.
I will take your word that $d_{iom}$ is normally distributed. Note that this is an assumption that often helps with the math, but is not necessarily representative. For example, you will always have some nonzero probability of a negative delay. However, if your expected value $d_i$ is significantly larger than $v_i$ that phenomenon can safely be neglected.
With that out of the way, I think you have at least 3 directions in which you could go, two scenario based, and chance-constrained optimization.
1. Monte Carlo Sampling
As Penghui Guo noted, you have all the information you need to sample possible realizations of $d_{iom}$: You assume it to be normally distributed, with center (i.e., expectation) $d_i$ and variance $v_i$. Now you can use your favorite random number generator to give you several samples, in C++ this could be std::normal_distribution. Each of the samples give you one scenario, you can weigh each scenario equally, but you will need many samples to describe the normal distribution appropriately.
2. Scenario-based modelling in 'plain-vanilla' stochastic programming style
Here you generate weighted scenarios, e.g., you know that 68.2% of all realizations are within $[d_i - 1v_i, d_i + 1v_i]$. You could represent this scenario with $d_{iom}=d_i$ to stay within your notation. Then you have one scenario with probability 15.9% to represent cases where $d_{iom} \in [-\inf, d_i - 1v_i]$, which you could represent with $d_{iom}=d_i-2v_i$, or $d_{iom}=d_i-1.5v_i$ (of course you can also compute the center of that tail of the distribution, which I did not bother to do). The third scenario would then represent the upper tail, so those cases where $d_{iom} \in [d_i + 1v_i, \inf]$, again with 15.9% probability.
A general note regarding scenarios
Scenarios blow up the model size incredibly fast! In particular as you have the three indices $i$, $o$, and $m$, which depending on your model, might all be dimensions in which you can generate scenarios. As a simple example, consider this: Let's say, you have $i \in {1,2}$ and $o \in {1,2}$ and we omit $m$. Now you already have 4 variables, $d_{1,1}$, $d_{1,2}$, $d_{2,1}$, $d_{2,2}$. If you only have 2 scenarios for each of these variables, you will end up with a scenario tree of 16 (=$2^4$) leafs. If you have 3 scenarios per variable, you'd already be at $3^4=81$ leafs! Thus, maybe the next framework is computationally more tractable and still an acceptable representation of reality:
3. Chance-constraint optimization
Let $\Phi(x)$ be the cumulative distribution function (or short cdf) of your delay, that means $\Phi(x) = \mathbb{P}(delay \leq x)$. This cdf depends on the variance $v_i$, so you can use $v_i$ to compute a value $d^{95}$ for $d_{iom}$ that covers, say, 95% of all delays. The exat percentage is your choice of course. This way your model is conservative in 95% (or 80%, or 98%,... you get the idea) of all cases, in the sense that it overestimates the delay and only in the remaining few cases you will have modeled the delay to be too short. Here is how you would model it:
$$d^{95} = \Phi^{-1}(95\%)$$
$$ c_i \geq s_i + d^{95}z_{iom}$$
where you'd only have to look up how to compute $\Phi^{-1}$ as a function of your variance $v_i$ in your computational environment. Essentially you are replacing your original inequality $ c_i \geq s_i + d_{iom}z_{iom}$ with the inequality $ \mathbb{P}(c_i \geq s_i + d_{iom}z_{iom}) \geq 95\%$
The charm of this solution is that you would only need one additional equation per delay for which you want chance constraints.
Hope this gives you some inspiration!
