Responding to Updated question:
The OP is properly creating a white noise process that is proportional to bandwidth, with the standard deviation increasing at $\sqrt{B}$ consistent with thermal noise given as $kTB$ where $k$ is Boltzmann's constant, $T$ is the temperature in Kelvin and $B$ is the bandwidth. We will assume for practical application that the thermal noise has been amplified to be sufficiently above the noise floor of the A/D converter (as the digital samples of the thermal noise).
I think the issue is in the way the OP is measuring the thermal noise. A proper measurement is the standard deviation of the signal, which (assuming the signal has zero mean, otherwise subtract the mean first) is simply the square root of the sum of the square of each sample. Doing this will result in the same standard deviation regardless of sampling rate as long as the bandwidth of the noise is higher than the sampling rate (white noise has infinite bandwidth but in practical applications it would be limited by the analog bandwidth of the A/D converter, or an anti-aliasing filter placed prior to A/D conversion).
So in summary, to measure the actual noise level as an rms quantity, take the average of all the samples and subtract this average from each sample (to be zero-mean). Then square each sample, sum the squares, and then take the square root. This is "root-mean-square" and for a zero mean signal will be the standard deviation of the signal. The square of the standard deviation is the variance. If the signal is in volts, then the variance divided by the resistance will be the power.
The OP was also interested (in the comments) in how to determine the maximum number of samples to use. For a white noise process, this is an excellent application of the Allan Deviation, available in python as the "allantools" package, with key points summarized in the plot below (plot is from this link with further details on the Allan Deviation for those interested):
Even if you don't have time to understand the full details of the Allan Deviation you can use the tool to create a plot as above for any random process using "frequency" as the input variable and not "phase". (The tool is intended for frequency variation vs time, and specifically a "White FM" process would be the same as any white noise process. Phase is an integral of frequency, so the tool allows for that optional input, but that would result in integrating the noise process which acts like a low pass filter-- just use "frequency" as the input and it will work as expected). The Allan Deviation shows a measure of standard deviation versus averaging time (it is more detailed than that as I explain in other posts but for white noise processes it is that simple). For stationary white noise specifically, the Allan Deviation will go down as the square root of the averaging time (consistent with the standard deviation of samples of white noise reducing by the square root on the number of samples averaged). For practical application many processes are only stationary (constant mean and standard deviation) over a short time interval, after which drift and other factors take over-- any computation with averaging time longer than this cross over point will result in a degraded measure of standard deviation. This has wide applicability to many processes where we are not time limited and want to know how long we should average for!
Details from prior answer below
When you double the sampling rate, you double the number of samples, which for a given threshold and same distribution, would double the number of samples crossing that threshold.
For example, below is a plot of a white noise process. To be "white noise" every sample would be independent of the next, so no matter how much we zoom in, this plot would look the same. If we set a threshold and sampled at twice the rate, over the same amount of time we would have twice the samples cross that threshold.
Thermal noise as we would sample with an A/D converter (ADC), if amplified to be higher than the ADC's own noise, is better represented as "band-limited white noise". This is because we would either low pass filter the signal before sampling (to eliminate all aliases), or if not, the ADC itself has a finite bandwidth. For this, as we zoom in, we will see the temporal variations in the waveform. Still if we sample this at twice the rate and have a threshold, we would again see twice the number of samples cross that threshold in the same amount of time.
I mention the above since thermal noise is a spectral density, meaning power over unit bandwidth in Hz. If the analog bandwidth increases, then the total power (and thus the standard deviation of the Gaussian distribution) would increase. This has nothing to do with the sampling rate but everything to do with the analog bandwidth that exists either from our own filtering or the finite bandwidth of the sampling device (real world).
For thermal noise, this would make sense only if the signal bandwidth prior to the analog to digital converter (thermal noise is an analog process) was higher than the sampling rate, and the thermal noise was amplified to be well above the quantization noise (assuming you were modelling what would be at the output of an A/D converter, otherwise we would be modelling the quantization noise from the A/D converter instead of the thermal noise: that would be very similar as a white noise process but with a uniform distribution instead of Gaussian).
Thermal noise at a fixed temperature is at a fixed power spectral density (given by $kT$ where $k$ is Boltzmann's constant and $T$ is the temperature in Kelvin) meaning the power over a unit bandwidth is constant. It is a white noise random process with a Gaussian (normal) distribution. At a given sampling rate, the total variance in time will equal the total variance in frequency (see Parseval's Theorem). This means if you create samples of a noise process using a Gaussian random generator with standard deviation $\sigma$ at sampling rate $f_s$, the power of that waveform ($\sigma^2)$ will be spread evenly over the frequency range given by the sampling rate (so the power spectral density will be $\sigma^2/f_s$). If the input to the ADC is properly low pass filtered such that there is no aliasing, then we would expect the total power due to amplified thermal noise to increase with increasing sampling rate (assuming we increase the bandwidth of the associated anti-alias filter to allow all the first Nyquist zone to pass). To model thermal noise at a constant power level in that case, you need to increase the $\sigma$ used by the sampling rate to account for the increased bandwidth. If we sample far below the analog bandwidth, we will get the same distribution (standard deviation and variance) independent of sampling rate, which can be explained by an aliasing model.
To help visualize this I have provided the graphic below. On the left are samples of a normalized ($\sigma = 1$) Gaussian white noise process. It is white because each sample is independent of all others (no memory) and when the mean is 0 as in this case, $\sigma$ is the root-mean-square or rms of the signal. On the right is a histogram, turned sideways to align with the amplitude distribution of the signal, and the histogram is showing the distribution of all amplitudes. A threshold is shown in red at $\sigma =1$. The point is, this waveform with a given $\sigma$ is independent of the sampling rate; I did not need to provide the sampling rate to define how this will look as long as each sample was independently generated. At any given sampling rate, with enough samples we will see the same distribution, and therefore over enough samples to be statistically valid will get the same percentage of samples crossing the threshold. If the time duration is the same, and we double the sampling rate, we will double the number of samples crossing the threshold. As the resulting digital samples, this is the same result we would get with noise filtered to the Nyquist bandwidth or noise with much higher bandwidth as long as the total power in both cases was the same (and we can see visually how we would get the same standard deviation for both those cases since it is based on the total power alone). If we were to provide filtering that is below the sampling rate, we would get different standard deviations depending on the bandwidth of the filter used. (This would be referred to as band-limited white noise).
