I am familiar and comfortable with the following definition of an integral (defined as the the Riemann Sum where n approaches infinity).
I have circled in red two terms to emphasize a concept that I was taught: the $\frac{b-a}{n}$ term effectively becomes the $dx$ term. This intuitively makes sense to me.
Additionally, I am familiar with the following statement (as expressed here: How does $u$-substitution work?)
However, I have boxed in red something that I find rather confusing when trying to interpret it through the understanding that $\frac{b-a}{n}$ sorta kinda equals $dx$...
Specifically, the $du$ value appears to be changing throughout the interval $[ g(a), g(b)]$. To reframe this in the "Riemann Sum" version,we would have something like: $\text{ "value that varies depending on where you are"} * \frac{g(b)-g(a)}{n}$
This is very confusing to me...because, if this is the correct interpretation, I'm not quite certain I understand why this is allowed.


You seem to be under the impression that one has to have one fixed $\Delta x$, i.e. a uniform subdivision of the interval. This is a very common misconception and the full definition of what a Riemann sum can be is given in the answer by LutzL.
If you have a Riemann sum with a uniform spacing then a u-substitution (performed at the level of the sum) will lead to a non-uniform spacing so you will have a $\Delta x$ that varies across the interval. For example consider $\int_0^1 2xe^{-x^2}{\rm d}x$ and let's use a uniform division of the interval $x_i = \frac{i}{n}$ to get $$R_n = \sum_{i=0}^{n-1} 2x_i e^{-x_i^2} \Delta x_i$$ where $\Delta x_i = x_{i+1}-x_i \equiv \frac{1}{n}$ is the same across the interval. In a u-substitution we will put $u_i = x_i^2$ to get $$R_n = \sum_{i=0}^{n-1} e^{-u_i} \Delta u_i \cdot \frac{x_i+x_1}{x_i+x_{i+1}}$$ where now $\Delta u_i = u_{i+1} - u_i = \frac{2i+1}{n^2}$ which varies from point to point. The last factor above will approach unity as $n\to\infty$ so for large $n$ this sum is just $\simeq \sum_{i=0}^{n-1} e^{-u_i} \Delta u_i$ which is a Riemann sum for the integral $\int_0^1 e^{-u}{\rm d}u$.
The largest spacing is found close to $u=1$ where $\Delta u_{n-1} = \frac{2n-1}{n^2} \sim \frac{2}{n}$ compared to the spacing close to $u=0$ where $\Delta u_0 = \frac{1}{n^2}$. However as $n\to\infty$ note that this largest spacing does go to zero and this is all we need for the Riemann sum to converge to the value of the integral.
Perhaps it's more instructive to show a figure. Recall that a Riemann sum (before taking the limit) is nothing but a finite sum of rectangles that approximate the (signed) area below the curve. Here is one example of such a sum where the spacing varies from point to point:
If we consider a finer subdivision of the interval, but still non-uniform, we generally get a much better approximation:
You can from this convince yourself that for approximating the area well it doesn't matter that some regions have a larger spacing than others as long as the largest spacing is "small enough". This is the basis of the general definition of a Riemann sum.